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ABSTRACT 


A  laboratory  spacecraft  simulator  testbed  is  first  introduced  to  examine  the 
problem  of  multiple  spacecraft  interacting  in  close  proximity.  This  testbed  enables 
validation  of  guidance,  navigation  and  control  (GNC)  algorithms  by  combining  6- 
Degrees  of  Freedom  (DoF)  computer  simulation  with  3-DoF  Hardware-In-the-Loop 
(HIL)  experimentation.  The  presented  3-DoF  spacecraft  simulator  employs  a  novel 
control  actuator  configuration  consisting  of  a  Miniature  Single  Gimbaled  Control 
Moment  Gyroscope  (MSGCMG)  and  dual  on/off  cold-gas  in-plane  vectorable  thrusters. 
The  dual  vectorable  thruster  design  enables  simultaneous  translation  and  attitude  control 
allowing  it  to  act  both  in  conjunction  with  the  MSGCMG  as  well  as  provide  sole  actuator 
control  throughout  a  commanded  closed-path  maneuver.  Small-time  local  controllability 
(STLC)  of  this  uniquely  actuated  system  via  Lie  Algebra  methods  is  formally 
demonstrated  and  results  of  experiments  conducted  on  the  described  testbed  are  included. 
From  this  study  in  3-DoF,  a  6-DoF  minimally  control  actuated  asymmetric  spacecraft 
design  is  proposed.  Six-DoF  control  of  this  underactuated  mechanical  system  is  achieved 
via  two  oppositely  mounted  hemispherically  vectorable  thrusters.  In  order  to  capitalize 
on  the  unique  nature  of  this  system  with  only  two  control  torques,  a  quaternion  feedback 
regulator  is  developed  to  yield  three-axis  stabilization  of  its  attitude.  This  regulator 
capitalizes  on  recent  advancements  in  generalized  inversion  and  perturbed  feedback 
linearizing  control  to  stabilize  the  dynamics  of  an  underactuated  asymmetric  spacecraft 
and  extends  this  to  include  stabilization  of  the  kinematics  of  the  system.  Two  control 
design  methodologies  are  derived.  The  first  is  Lyapunov  based,  yielding  a  globally  stable 
system,  while  the  second  yields  local  stability  within  a  domain  of  attraction  through 
perturbed  feedback  linearization.  Results  of  several  numerical  simulations  are  presented 
for  an  asymmetric  spacecraft  with  two  bounded  body-fixed  control  torques.  The 
proposed  attitude  control  method  is  not  intended  to  provide  attitude  maintenance  for 
attitude  tracking  or  in  the  presence  of  relatively  large  disturbance  torques;  however,  it 
may  prove  widely  applicable  for  detumbling  and  reorientation  maneuvers  of  spacecraft 
with  only  two  available  control  torques. 
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I. 


INTRODUCTION 


The  traditional  spacecraft  system  is  a  monolithic  structure  with  a  single  mission 
focused  design  and  lengthy  production  and  qualification  schedules  coupled  with 
enormous  cost.  Additionally,  there  rarely,  if  ever,  are  any  designed  preventive 
maintenance  plans  or  re-fueling  capability.  There  has  been  much  research  in  recent  years 
into  alternative  options.  One  option  involves  autonomous  on-orbit  servicing  of  current  or 
future  monolithic  spacecraft  systems.  The  U.S.  Department  of  Defense  (DoD)  embarked 
on  a  highly  successful  venture  to  prove  such  a  concept  with  the  Defense  Advanced 
Research  Projects  Agency’s  (DARPA)  Orbital  Express  program.  Orbital  Express 
demonstrated  all  of  the  enabling  technologies  required  for  autonomous  on-orbit  servicing 
to  include  refueling,  component  transfer,  autonomous  satellite  grappling  and  berthing, 
rendezvous,  inspection,  proximity  operations,  docking  and  undocking,  and  autonomous 
fault  recognition  and  anomaly  handling  (Kennedy  2008).  Another  potential  option 
involves  a  paradigm  shift  from  the  monolithic  spacecraft  system  to  one  involving 
multiple  interacting  spacecraft  that  can  autonomously  assemble  and  reconfigure. 
Numerous  benefits  are  associated  with  autonomous  spacecraft  assemblies,  ranging  from  a 
removal  of  significant  intra-modular  reliance  that  provides  for  parallel  design, 
fabrication,  assembly  and  validation  processes  to  the  inherent  smaller  nature  of 
fractionated  systems  that  allows  for  each  module  to  be  placed  into  orbit  separately  on 
more  affordable  launch  platforms  (Mathieu  and  Weigel  2005). 

A.  GROUND-BASED  HARDWARE-IN-THE-LOOP  EXPERIMENTAL 

VALIDATION 

With  respect  specifically  to  the  validation  process,  the  significantly  reduced 
dimensions  and  mass  of  aggregated  spacecraft  when  compared  to  the  traditional 
monolithic  spacecraft  allow  for  not  only  component  but  even  full-scale  on-the-ground 
Hardware-In-the-Loop  (HIL)  experimentation.  Likewise,  much  of  the  HIE 
experimentation  required  for  on-orbit  servicing  of  traditional  spacecraft  systems  can  also 
be  accomplished  in  ground-based  laboratories  (Creamer  2007).  This  type  of  HIL 
experimentation  complements  analytical  methods  and  numerical  simulations  by 
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providing  a  low-risk,  relatively  low-cost  and  potentially  high-return  method  for 
validating  the  technology,  navigation  techniques  and  control  approaches  associated  with 
spacecraft  systems.  Several  approaches  exist  for  the  actual  HIL  testing  in  a  laboratory 
environment  with  respect  to  spacecraft  guidance,  navigation  and  control.  One  such 
method  involves  reproduction  of  the  kinematics  and  vehicle  dynamics  for  3 -DoF  (two 
horizontal  translational  degrees  and  one  rotational  degree  about  the  vertical  axis)  through 
the  use  of  robotic  spacecraft  simulators  that  float  via  planar  air  bearings  on  a  flat 
horizontal  floor.  This  particular  method  is  currently  being  employed  by  several  research 
institutions  and  is  the  validation  method  of  choice  for  our  research  into  GNC  algorithms 
for  proximity  operations  at  the  Naval  Postgraduate  School  (Machida,  Toda,  and  Iwata 
1992;  Ullman  1993;  Corrazzini  and  How  1998;  Marchesi,  Angrilli  and  Venezia  2000; 
Ledebuhr  et  al.  2001;  Nolet,  Kong,  and  Miller  2005;  LeMaster,  Schaechter,  and 
Carrington  2006;  Romano,  Friedman,  and  Shay  2007). 

For  spacecraft  involved  in  proximity  operations,  the  in-plane  and  cross-track 
dynamics  are  decoupled,  as  modeled  by  the  Hill-Clohessy-Wiltshire  (HCW)  equations, 
thus  the  reduction  to  3-Degree  of  Freedom  (DoF)  does  not  appear  to  be  a  critical  limiter. 
One  consideration  involves  the  reduction  of  the  vehicle  dynamics  to  one  of  a  double 
integrator.  However,  the  orbital  dynamics  can  be  considered  to  be  a  disturbance  that 
needs  to  be  compensated  for  by  the  spacecraft  navigation  and  control  system  during  the 
proximity  navigation  and  assembly  phase  of  multiple  systems.  Thus,  the  flat  floor 
testbed  can  be  used  to  capture  many  of  the  critical  aspects  of  an  actual  autonomous 
proximity  maneuver  that  can  then  be  used  for  validation  of  numerical  simulations. 
Portions  of  the  here-in  described  testbed,  combined  with  the  first  generation  robotic 
spacecraft  simulator  of  the  Spacecraft  Robotics  Laboratory  (SRL)  at  the  Naval 
Postgraduate  School  (NPS),  have  been  employed  to  propose  and  experimentally  validate 
control  algorithms.  The  interested  reader  is  referred  to  Romano  et  al.  (2007)  for  a  full 
description  of  this  robotic  spacecraft  simulator  and  the  associated  HIL  experiments 
involving  its  demonstration  of  successful  autonomous  spacecraft  approach  and  docking 
maneuvers  to  a  collaborative  target  with  a  prototype  docking  interface  of  the  Orbital 
Express  program. 
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A  unique  control  problem  exists,  given  a  requirement  for  spacecraft  aggregates  to 
rendezvous  and  dock  during  the  final  phases  of  assembly  and  a  desire  to  maximize  the 
useable  surface  area  of  the  spacecraft  for  power  generation,  sensor  packages,  docking 
mechanisms  and  payloads  while  minimizing  thruster  impingement.  In  fact,  control  of 
such  systems  using  the  standard  control  actuator  configuration  of  fixed  thrusters  on  each 
face  coupled  with  momentum  exchange  devices  can  be  challenging,  if  not  impossible. 
For  such  systems,  a  new  and  unique  configuration  is  proposed  consisting  of  vectorable 
thrusters  that  may  capitalize,  for  instance,  on  the  recently  developed  carpal  robotic  joint 
invented  by  Canfield  and  Reinholtz  (1998)  with  its  hemispherical  vector  space.  It  is  here 
demonstrated  through  Lie  algebra  analytical  methods  and  experimental  results  that  two 
vectorable  in-plane  thrusters  in  an  opposing  configuration  can  yield  a  minimum  set  of 
actuators  for  a  controllable  reduced  order  (3-DoF)  system.  It  will  be  shown  that  by 
coupling  the  proposed  set  of  vectorable  in-plane  thrusters  with  a  single  degree  of  freedom 
torquer  such  as  a  Control  Moment  Gyroscope,  an  additional  degree  of  redundancy  can  be 
gained  for  the  reduced  order  system.  Experimental  results  are  included  using  SRL’s 
second-generation,  reduced-order  spacecraft  simulator  with  a  state  feedback  linearized 
controller  to  demonstrate  its  ability  to  navigate  a  closed  circular  path  with  the  proposed 
actuator  configuration.  A  general  overview  of  this  spacecraft  simulator  is  presented  with 
additional  details  on  the  simulators  being  found  in:  Hall  (2006),  Eikenberry  (2006),  Price 
(2006),  Romano  and  Hall  (2006),  Hall  and  Romano  (2007a),  and  Hall  and  Romano 
(2007b). 

B.  ROT O-TRAN SL ATION  OF  AN  UNDERACTUATED  SPACECRAFT 

With  respect  to  the  full  order  (6-DoF)  system  of  roto-translation,  it  will  be 
demonstrated  that  a  set  of  two  vectorable  hemispherical  thrusters  in  an  opposing 
configuration  can  yield  a  minimum  set  of  actuators  for  a  controllable  relative  motion 
spacecraft  system.  The  proposed  actuator  configuration  can  readily  be  seen  to  yield  an 
underactuated  mechanical  system,  that  is  to  say  the  number  of  controls  Nu  is  fewer  than 

the  degrees  of  freedom.  It  will  be  shown  that  by  adding  a  single-degree  of  freedom 
torque,  the  system  can  become  fully  actuated  and  a  state-feedback  linearizing  controller 
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can  be  designed  to  drive  the  system  to  the  desired  position  and  orientation.  Without  the 
single-degree  of  freedom  torque,  this  control  actuator  configuration  equates  to  only  two 
of  the  three  body-fixed  axes  being  directly  torque  controlled  while  the  third  can  be 
indirectly  controlled  by  capitalizing  on  the  coupling  of  the  terms  of  the  moments  of 
inertia  for  the  system  that  appear  in  the  Euler  equations.  Furthermore,  given  the  cascaded 
nature  of  the  angular  motion  equations  whereby  the  angular  rates  drive  the  orientation, 
there  has  been  much  research  to  date  to  develop  control  algorithms  to  provide  three-axis 
stabilization  for  an  underactuated  spacecraft  with  only  two  control  torques.  However,  to 
date,  there  has  not  been  a  smooth  time-invariant  control  algorithm  to  provide  attitude 
stabilization  of  such  a  spacecraft  with  arbitrary  inertia.  By  considering  the  attitude 
stabilization  for  the  underactuated  spacecraft  system  of  angular  motion  equations  in  the 
general  case,  it  may  prove  to  be  widely  applicable  to  not  only  proximity  operations  but 
also  de-tumbling  and  reorientation  maneuvers  of  underactuated  spacecraft  that  may  either 
be  designed  with  only  two  directly  actuated  control  axes  or  be  experiencing  control 
actuator  failure  about  one  of  its  control  axes  during  their  mission  life. 

C.  STATE-OF-THE-ART  IN  UNDERACTUATED  RIGID  BODY 

STABILIZATION 

The  problem  of  stabilization  of  a  rigid  spacecraft’s  attitude  dynamics  and 
kinematics  has  been  studied  over  the  years  in  many  papers  and  articles.  However,  the 
vast  majority  of  the  proposed  control  laws  assume  that  the  spacecraft  is  fully  actuated. 
Wie  and  Barba  (1985),  Wie,  Weiss  and  Arapostathis  (1989),  Vadali  (1989),  and  Bajodah 
(2009a)  address  several  nonlinear  control  techniques  that  provide  time-invariant  global 
asymptotic  stability  of  the  fully  actuated  spacecraft  system  of  equations.  Although  these 
control  laws  provide  for  the  necessary  control  of  a  nominally  designed  three-axis 
stabilized  spacecraft  in  which  three  control  torques  exist,  the  question  of  control  of 
underactuated  spacecraft  naturally  enters  when  discussing  actuator  failures  or  when 
proposing  minimally  designed  spacecraft  systems.  It  is  well  understood  that  full  order  (3- 
DoF)  control  of  the  kinematics  of  such  underactuated  systems  presents  a  challenging 
control  problem;  however,  it  should  also  be  recognized  that  it  has  the  distinct  potential  to 
provide  several  key  benefits.  Specifically,  under  the  present  thrust  of  Operationally 


4 


Responsive  Space,  where  one  focus  is  on  transitioning  from  the  typical  large  monolithic 
spacecraft  design  to  one  that  embraces  the  new  spacecraft  paradigm  of  smaller,  faster  to 
manufacture  and  cheaper  to  produce  and  employ,  underactuated  control  could  provide  a 
key  enabling  technology.  Furthermore,  in  light  of  many  current  traditional  spacecraft 
systems  remaining  in  operation  well  past  their  intended  design  life  despite  actuator 
failures  that  degrade  their  capabilities,  underactuated  control  could  enable  these  aging 
systems  to  satisfy  their  original  missions. 

The  investigation  of  stabilization  of  underactuated  spacecraft  kinematics  and 
attitude  dynamics  began  with  the  theoretical  establishment  of  the  necessary  and  sufficient 
conditions  for  the  controllability  of  a  rigid  body’s  attitude  with  either  gas  thrusters  or 
momentum  exchange  devices  by  Crouch  (1984).  He  concluded  that,  for  a  spacecraft  with 
momentum  exchange  devices,  controllability  is  impossible  with  fewer  than  three,  while 
for  a  spacecraft  with  independent  paired  jets,  controllability  is  possible  with  two.  It  was 
later  demonstrated  by  Kerai  (1995),  by  using  geometric  control  theory,  that  small-time 
local  controllability  of  the  rigid  body  equations  assuming  paired  gas  jets  can  indeed  be 
achieved  with  only  two  control  torques.  Byrnes  and  Isidori  (1991)  proved  that  the  full 
angular  motion  equations  for  a  rigid  spacecraft  with  only  two  controls  cannot  be 
asymptotically  stabilized  by  smooth  pure  state  feedback  because  they  violate  Brockett’s 
(1983)  theorem  on  non-holonomic  underactuated  systems.  With  this  in  mind,  they 
proposed  a  smooth  feedback  controller  to  affect  partial  stabilization  of  the  rigid  body 
model  resulting  in  a  revolute  constant-rate  motion  about  the  uncontrolled  axis  of  rotation. 
Later,  Krishnan,  McClamroch  and  Reyhanoglu  (1995)  proposed  a  hybrid  control  design 
combining  continuous  time  features  with  discrete  event  features  to  affect  a  discontinuous 
feedback  control  strategy  to  stabilize  any  equilibrium  attitude  of  an  underactuated 
spacecraft  with  two  momentum  wheel  actuators  in  finite  time  under  the  restriction  that 
the  total  angular  momentum  vector  of  the  system  is  zero.  This  control  methodology 
translates  directly  to  a  study  of  an  underactuated  axi-symmetric  spacecraft.  Tsiotras, 
Corless,  and  Longuski  (1995)  and  Tsiotras  and  Luo  (2000)  also  dealt  with  control  of 
underactuated  axi-symmetric  spacecraft  by  proposing  a  time-invariant  feedback  control 
law  to  asymptotically  stabilize  the  orientation  of  two  of  the  three  body-fixed  axes.  In 
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addition  to  only  providing  for  partial  attitude  stabilization  of  axi-symmetric  spacecraft, 
their  discussion  was  limited  to  cases  where  the  angular  velocity  about  the  unactuated  axis 
is  zero  at  the  start  of  the  maneuver.  Tsiotras  and  Schleicher  (2000),  and  Tsiotras  and 
Doumtchenko  (2000)  relaxed  the  restriction  on  the  symmetry  of  the  spacecraft  slightly  to 
consider  a  nearly  axi-symmetric  spacecraft  by  a  small  parameter  and  a  set  of  time- 
invariant  control  laws  are  proposed  to  stabilize  the  angular  velocity  and  attitude  of  a 
spacecraft  about  a  certain  axis  by  virtual  control  inputs  of  the  two  actuated  angular  rates. 

The  global  asymptotic  rate  stabilization  problem  without  concern  for  kinematics 
of  a  fully  asymmetric  underactuated  rigid  spacecraft  was  addressed  by  Coverstone- 
Carroll  (1996)  through  the  use  of  a  Variable  Structure  Controller  (VSC).  Bajodah 
(2009b)  also  addressed  the  rate-only  stabilization  problem  for  detumbling  maneuvers 
through  the  use  of  singularly  perturbed  feedback  linearization  and  generalized  inverse 
control  methodologies.  Although  both  of  these  controllers  prove  to  be  robust  to  large 
initial  angular  velocities  around  all  three  axes  in  the  presence  of  actuator  torque 
limitations,  they  both  require  an  additional  controller  to  provide  desired  kinematic 
alignment  after  the  detumbling  maneuver.  One  such  controller,  as  proposed  by 
Coverstone-Carroll  (1996),  is  a  simple  linear  controller  that  is  used  to  perform  a  series  of 
eigenaxis  rotations  which  precludes  smooth  attitude  tracking. 

The  problem  of  stabilization  of  both  the  kinematics  and  dynamics  of  an 
underactuated  asymmetric  spacecraft  was  most  recently  addressed  by  Casagrande, 
Astolfi,  and  Parisini  (2008)  who  proposed  a  time -variant  switching  control  law  to  effect 
global  asymptotic  stabilization  of  the  closed-loop  system.  Although  novel,  the  proposed 
law  lacks  detailed  simulation  results  by  considering  only  the  case  where  the  initial 
angular  rates  about  two  of  the  axes  to  include  the  unactuated  axis  are  initially  zero. 
Furthermore,  real-world  spacecraft  with  flexible  parts,  antennas,  fuel  slosh,  etc.,  may 
preclude  the  use  of  time -variant  control  laws  because  they  have  the  distinct  potential  of 
producing  unacceptable  transient  response  and  might  therefore  lead  to  instability 
(Tsiotras  and  Doumtchenko,  2000).  Behai  et  al.  (2002)  address  the  nonlinear  tracking 
control  of  an  axi-symmetric  spacecraft  by  developing  a  kinematic  controller  to  determine 
the  desired  actuated  angular  rates  which  are  in  turn  used  as  control  inputs  to  the  dynamic 
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system  through  the  use  of  standard  back-stepping  techniques.  This  method  yields  only 
asymptotic  dynamic  and  kinematic  stabilization  results  for  an  axi-symmetric  rigid  body 
given  the  restriction  that  the  angular  rate  about  the  unactuated  axis  is  close  to  zero  but  it 
does  yield  bounded  results  otherwise. 

The  goal  of  this  work  is  to  extend  the  research  into  control  algorithms  for 
underactuated  rigid  spacecraft  attitude  control  by  proposing  a  novel,  time-invariant 
smooth  quaternion  feedback  regulator  based  on  generalized  inverse  methods  to  affect 
three-axis  attitude  stabilization  of  the  error  quaternion  kinematics  for  an  underactuated 
rigid  spacecraft  with  arbitrary  inertia  matrices  and  two  realistically  bounded  body-fixed 
torques  for  required  reorientation  maneuvers.  The  problem  of  three-axis  stabilization  of 
the  attitude  of  an  underactuated  spacecraft  with  direct  control  about  only  two  of  the  body- 
fixed  control  axes  will  be  addressed  in  the  general  case  where  the  spacecraft’s  attitude  is 
referenced  to  an  inertially  fixed  frame.  From  this,  the  proposed  quaternion  feedback 
regulator  can  be  shown  to  be  seamlessly  modified  to  account  for  attitude  control  with 
respect  to  a  relative  motion  frame  of  reference  as  is  the  case  for  a  chaser-target 
rendezvous  maneuver.  After  affecting  attitude  error  stabilization,  a  spacecraft  can  be 
propelled  towards  another  spacecraft  via  various  navigation  schemes  such  as 
conventional  waypoint  navigation. 

D.  SCOPE  OF  THE  DISSERTATION 

This  dissertation  advances  the  body  of  knowledge  with  respect  to  control  of 
underactuated  spacecraft  in  three  key  areas: 

1.  Laboratory  experimentation  of  a  reduced-order  underactuated  spacecraft 
simulator  with  vectorable  thrusters  and  a  miniature  control  moment 
gyroscope.  Using  feedback  linearizing  control  methodology  coupled  with 
Schmitt  Trigger  and  Pulse  Width  Modulation  logic,  experimental  results  are 
presented  which  validate  the  capability  of  this  novel  control  actuator  design  to 
propel  the  spacecraft  simulator  around  a  tightly  constrained  path. 

2.  Analytical  detennination  of  the  small-time  local  controllability  of  a  generic 
full-order  spacecraft  under  variations  on  the  control  inputs.  This  study  is  able 


7 


to  provide  the  interested  spacecraft  systems  engineer  with  the  ability  to 
determine  the  minimum  number  of  control  actuators  necessary  to  maintain 
controllability.  Furthermore,  this  analysis  can  be  used  to  aid  in  dealing  with 
both  control  actuator  failures  on  existing  spacecraft  systems  or  planning  for 
minimally  designed  spacecraft. 

3.  A  smooth  time-invariant  state  feedback  control  logic  based  on  quaternion 
feedback  regulation  is  derived  to  yield  stabilization  of  the  error  kinematics  of 
a  spacecraft  with  only  two  control  torques  and  arbitrary  inertia.  Two  separate 
control  designs  are  presented,  the  first  being  Lyapunov  function  based  and  the 
other  being  perturbed  feedback  linearizing  in  nature.  Results  of  the  numerical 
results  considering  both  of  these  designs  are  presented  for  various  maneuvers. 
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II.  LABORATORY  EXPERIMENTATION  OF  GUIDANCE  AND 
CONTROL  OF  SPACECRAFT  DURING  PROXIMITY 

MANEUVERS 

While  presenting  an  overview  of  a  robotic  testbed  for  HIL  experimentation  of 
guidance  and  control  algorithms  for  on-orbit  proximity  maneuvers,  this  chapter 
specifically  focuses  on  exploring  the  feasibility,  design  and  evaluation  in  a  3-DoF 
environment  of  a  vectorable  thruster  configuration  combined  with  optional  miniature 
single  gimbaled  control  moment  gyro  (MSGCMG)  for  an  agile  small  spacecraft. 
Specifically,  the  main  aims  are  to  present  and  practically  confirm  the  theoretical  basis  of 
small-time  local  controllability  for  this  unique  actuator  configuration  through  both 
analytical  and  numerical  simulations  performed  in  previous  works  (Romano  and  Hall 
2006;  Hall  and  Romano  2007a;  Hall  and  Romano  2007b)  and  to  validate  the  viability  of 
using  this  minimal  control  actuator  configuration  on  a  small  spacecraft  in  a  practical  way. 
Furthermore,  the  experimental  work  is  used  to  confirm  the  controllability  of  this 
configuration  along  a  fully  constrained  trajectory  through  the  employment  of  a  smooth 
feedback  controller  based  on  state  feedback  linearization  and  linear  quadratic  regulator 
techniques  and  proper  state  estimation  methods.  The  chapter  is  structured  as  follows: 
First  the  design  of  the  experimental  testbed  including  the  floating  surface  and  the  second 
generation  3-DoF  spacecraft  simulator  is  introduced.  Then  the  dynamics  model  for  the 
spacecraft  simulator  with  vectorable  thrusters  and  momentum  exchange  device  are 
formulated.  The  controllability  associated  with  this  uniquely  configured  system  is  then 
addressed  with  a  presentation  of  the  minimum  number  of  control  inputs  to  ensure  small 
time  local  controllability.  Next,  a  fonnal  development  is  presented  for  the  state  feedback 
linearized  controller,  state  estimation  methods,  Schmitt  trigger  and  Pulse  Width 
Modulation  scheme.  Finally,  experimental  results  are  presented  that  demonstrate  a 
closed-path  circular  trajectory  about  an  arbitrary  reference  that  is  representative  of  a 
possible  inspection  of  a  target  spacecraft  by  a  given  chaser. 
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A. 


THE  NPS  ROBOTIC  SPACECRAFT  SIMULATOR  TESTBED 


Three  generations  of  robotic  spacecraft  simulators  have  been  developed  at  the 
NPS  Spacecraft  Robotics  Laboratory,  in  order  to  provide  for  relatively  low-cost  HIL 
experimentation  of  GNC  algorithms  for  spacecraft  proximity  maneuvers  (see  Figure  1). 
In  particular,  the  second  generation  robotic  spacecraft  simulator  testbed  is  used  for  the 
here-in  presented  research.  The  whole  spacecraft  simulator  testbed  consists  of  three 
components.  The  two  components  specifically  dedicated  to  HIL  experimentation  in  3- 
DoF  are  a  floating  surface  with  an  indoor  pseudo-GPS  (iGPS)  measurement  system  and 
one  3 -DoF  autonomous  spacecraft  simulator.  The  third  component  of  the  spacecraft 
simulator  testbed  is  a  6-DoF  simulator  stand-alone  computer  based  spacecraft  simulator 
and  is  separated  from  the  HIL  components.  Additionally,  an  off-board  desktop  computer 
is  used  to  support  the  3-DoF  spacecraft  simulator  by  providing  the  capability  to  upload 
software,  initiate  experimental  testing,  receive  logged  data  during  testing  and  process  the 
iGPS  position  coordinates.  Figure  2  depicts  the  robotic  spacecraft  simulator  in  the 
Proximity  Operations  Simulator  Facility  (POSF)  at  NPS  with  key  components  identified. 
The  main  testbed  systems  are  briefly  described  in  the  next  sections  with  further  details 
given  in  Hall  (2006),  Price  (2006),  Eikenberry  (2006),  Romano  and  Hall  (2006),  Hall  and 
Romano  (2007a),  and  Hall  and  Romano  (2007b). 
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Figure  1.  Three  Generations  of  Spacecraft  Simulators  at  the  NPS  Spacecraft 
Robotics  Laboratory  (first,  second  and  third  generations,  from  left  to  right) 


1.  Floating  Surface 

A  4.9  m  by  4.3  m  epoxy  floor  surface  provides  the  base  for  the  floatation  of  the 
spacecraft  simulator.  The  use  of  planar  air  bearings  on  the  simulator  reduces  the  friction 

■j 

to  a  negligible  level  and  with  an  average  residual  slope  angle  of  approximately  2.6x10' 
deg  for  the  floating  surface,  the  average  residual  acceleration  due  to  gravity  is 
approximately  1.8x10'  ms'  .  This  value  of  acceleration  is  two  orders  of  magnitude  lower 
than  the  nominal  amplitude  of  the  measured  acceleration  differences  found  during 
reduced  gravity  phases  of  parabolic  flights  (Romano  et  al.  2007). 
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Figure  2.  SRL’s  2nd  Generation  3 -DoF  Spacecraft  Simulator 


2.  Three-DoF  Robotic  Spacecraft  Simulator 

SRL’s  second  generation  robotic  spacecraft  simulator  is  modularly  constructed 
with  three  easily  assembled  sections  dedicated  to  each  primary  subsystem.  Prefabricated 
6105-T5  Aluminum  fractional  t-slotted  extrusions  form  the  cage  of  the  vehicle  while  one 
square  foot,  0.25  inch  thick  static  dissipative  rigid  plastic  sheets  provide  the  upper  and 
lower  decks  of  each  module.  The  use  of  these  materials  for  the  basic  structural 
requirements  provides  a  high  strength  to  weight  ratio  and  enable  rapid  assembly  and 
reconfiguration.  Table  1  reports  the  key  parameters  of  the  3-DoF  spacecraft  simulator. 
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Subsystem 

Characteristic 

Parameter 

Structure 

Length  and  width 

Height 

Mass 

J, 

0.30  m 

0.69  m 

26  kg 

0.40  kg-m2 

Propulsion 

Propellant 

Equiv.  storage  capacity 
Operating  pressure 

Thrust  (x2) 

ISP 

Total  AV 

Compressed  Air 

0.05  m3  @31.03  MPa 

0.41  MPa 

0.159  N 

34.3  s 

31.1  m/s 

Flotation 

Propellant 

Equiv.  storage  capacity 
Operating  pressure 

Linear  air  bearing  (x4) 
Continuous  operation 

Air 

0.05  m3  @31.03  MPa 
.51  MPa 

32  mm  diameter 
-40  min 

CMG  Attitude  Control 

Max  torque 

Momentum  storage 

0.668  Nm 

0.098  Nms 

Electrical  &  Electronic 

Battery  type 

Storage  capacity 
Continuous  Operation 
Computer 

Lithium-Ion 

12  Ah  @  28V 
~6  h 

1  PC  104  Pentium  III 

Sensors 

Fiber  optic  gyro 

Position  sensor 
Magnetometer 

KVH  Model  DSP-3000 
Metris  iGPS 

MicroStrain  3DM-GX1 

Table  1.  Key  Parameters  of  the  2nd  Generation  3-DoF  Robotic  Spacecraft 

Simulator 


a.  Propulsion  and  Flotation  Subsystems 

The  lowest  module  houses  the  flotation  and  propulsion  subsystems.  The 
flotation  subsystem  is  composed  of  four  planar  air  bearings,  an  air  filter  assembly,  dual 
4500  PSI  (31.03  MPa)  carbon-fiber  spun  air  cylinders  and  a  dual  manifold  pressure 
reducer  to  provide  75  PSI  (0.5 1  MPa).  This  pressure  with  a  volume  flow  rate  for  each  air 
bearing  of  3.33  slfim  (3.33  x  10’  m  /min)  is  sufficient  to  keep  the  simulator  in  a  friction- 
free  state  for  nearly  40  minutes  of  continuous  experimentation  time.  The  propulsion 
subsystem  is  composed  of  dual  vectorable  supersonic  on-off  cold-gas  thrusters  and  a 
separate  dual  carbon-fiber  spun  air  cylinder  and  pressure  reducer  package  regulated  at  60 
PSI  (0.41  MPa)  and  has  the  capability  of  providing  the  system  31.1  m/s  AV . 
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b.  Electronic  and  Power  Distribution  Subsystems 

The  power  distribution  subsystem  is  composed  of  dual  lithium-ion 
batteries  wired  in  parallel  to  provide  28  volts  for  up  to  12  Amp-Hours  and  is  housed  in 
the  second  deck  of  the  simulator.  A  four  port  DC-DC  converter  distributes  the  requisite 
power  for  the  system  at  5,  12  or  24  volts  DC.  An  attached  cold  plate  provides  heat 
transfer  from  the  array  to  the  power  system  mounting  deck  in  the  upper  module.  The 
current  power  requirements  include  a  single  PC- 104  CPU  stack,  a  wireless  router,  three 
motor  controllers,  three  separate  normally-closed  solenoid  valves  for  thruster  and  air 
bearing  actuation,  a  fiber  optic  gyro,  a  magnetometer  and  a  wireless  server  for 
transmission  of  the  vehicle’s  position  via  the  pseudo-GPS  system. 

c.  Translation  and  Attitude  Control  System  Actuators 

The  3-DoF  robotic  spacecraft  simulator  includes  actuators  to  provide  both 
translational  control  and  attitude  control.  A  full  development  of  the  controllability  for 
this  unique  configuration  of  dual  rotating  thrusters  and  one-axis  Miniature-Single 
Gimbaled  Control  Moment  Gyro  (MSGCMG)  will  be  demonstrated  in  subsequent 
sections  of  this  paper.  The  translational  control  is  provided  by  two  cold-gas  on-off 
supersonic  nozzle  thrusters  in  a  dual  vectorable  configuration.  Each  thruster  is  limited  in 
a  region  ±n/2  with  respect  to  the  face  nonnal  and,  through  experimental  testing  at  the 
supplied  pressure,  has  been  demonstrated  to  have  an  ISP  of  34.3  s  and  able  to  provide 
0.159  N  of  thrust  with  less  than  10  msec  actuation  time  (Lugini  and  Romano  2009).  The 
MSGCMG  is  capable  of  providing  0.668  Nm  of  torque  with  a  maximum  angular 
momentum  of  0.098  Nms. 

3.  Six-DoF  Computer-based  Numerical  Spacecraft  Simulator 

A  separate  component  of  SRL’s  spacecraft  simulator  testbed  at  NPS  is  a  6-DoF 
computer-based  spacecraft  simulator.  This  simulator  enables  full  6-DoF  numerical 
simulations  to  be  conducted  with  realistic  orbital  perturbations  including  aerodynamic, 
solar  pressure  and  third-body  effects,  and  earth  oblateness  up  to  J4.  Similar  to  the  3-DoF 
robotic  simulator,  the  numerical  simulator  is  also  modularly  designed  within  a 
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MATLAB®/Simulink®  architecture  to  allow  near  seamless  integration  and  testing  of 
developed  guidance  and  control  algorithms.  Additionally,  by  using  the 
MATLAB®/Simulink®  architecture  with  the  added  Real  Time  Workshop™  toolbox,  the 
developed  control  algorithms  can  be  readily  transitioned  into  C-code  for  direct 
deployment  onto  the  3 -DoF  robotic  simulator’s  onboard  processor.  A  full  discussion  of 
the  process  by  which  this  is  accomplished  and  simplified  for  rapid  real-time 
experimentation  on  the  3 -DoF  testbed  for  either  the  proprietary  MATLAB®  based 
XPCTarget™  operating  system  is  given  in  (Hall  2006;  Price  2006)  or  for  an  open-source 
Linux  based  operating  system  with  the  Real  Time  Application  Interface  (RTAI)  is  given 
in  (Bevilacqua,  Hall,  Horning,  and  Romano  2009). 

B.  DYNAMICS  OF  A  3-DOF  SPACECRAFT  SIMULATOR  WITH 
VECTORABLE  THRUSTERS  AND  MOMENTUM  EXCHANGE  DEVICE 

Two  sets  of  coordinate  frames  are  established  for  reference:  an  inertial  frame  7 

with  orthogonal  axes  defined  by  the  unit  vectors  j  /, ,  i2 ,  /)  j  and  body-fixed  frame  7>  with 

orthogonal  axes  defined  by  the  unit  vectors  .  These  reference  frames  are 

depicted  in  Figure  4  along  with  the  necessary  external  forces  and  parameters  required  to 
properly  define  the  simulators  motion.  The  origin  of  the  body-fixed  coordinate  system  is 
taken  to  be  the  center  of  mass  C  of  the  spacecraft  simulator  and  is  assumed  to  be 

collocated  with  the  simulator’s  geometric  center.  By  the  effects  of  the  flat  floor,  6,  is 
maintained  aligned  with  i3  while  b{  is  defined  to  be  in  line  with  the  thrusters  points  of 
action.  The  position  and  velocity  vectors  of  2?  with  respect  to  7  expressed  in  2?  are 
given  by  r  and  v  so  that  rBI  marks  the  position  of  the  simulator  with  respect  to  the  origin 

of  7  as  measured  by  the  inertial  measurement  sensors  and  provides  the  vehicle’s  two 
degrees  of  translational  freedom.  The  vehicle’s  rotational  freedom  is  described  by  an 

angle  of  rotation  <9,  between  bx  and  /,  about  4  .  The  angular  velocity  of  7>  with  respect 
to  7  expressed  in  7  is  thus  limited  to  one  degree  of  freedom  and  is  denoted  by  co3 .  The 
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spacecraft  simulator  is  assumed  to  be  rigid  and  therefore  a  constant  moment  of  inertia 
( )  exists  about  the  6, .  Furthennore,  any  changes  to  the  mass  of  the  simulator  ( m  )  due 
to  thruster  firing  are  neglected. 

The  forces  imparted  at  a  distance  L  from  the  center  of  mass  by  the  vectorable  on- 
off  thrusters  are  denoted  by  Fa  and  Fh  respectively.  The  thrust  angle  au  defines  the 

orientation  of  thruster  a  in  22  and  is  the  angle  measured  from  b]  in  a  clockwise 
direction  (right-hand  rotation)  to  Ffl .  Likewise,  thrust  angle  ab  defines  the  orientation  of 
thruster  b  in  2?  and  is  the  angle  measured  from  -bx  in  a  clockwise  direction  (right-hand 
rotation)  to  .  The  torque  imparted  on  the  vehicle  by  a  momentum  exchange  device 
such  as  the  MSGCMG  is  denoted  by  Tcm  and  can  be  constrained  to  exist  only  about  the 
b3  axis. 


Figure  3.  Schematic  of  SRL’s  2nd  Generation  Spacecraft  Simulator 


1,  Controllability  Analysis  of  a  Single-Gimbaled  Control  Moment 
Gyroscope  for  a  Reduced  Order  System 

The  feasibility  and  utility  of  using  a  single-gimbaled  control  moment  gyroscope 
for  the  actuation  of  a  3-DoF  spacecraft  simulator  can  easily  be  demonstrated  through 
development  of  the  reduced  order  angular  motion  equations  for  the  vehicle  and  the 
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MSGCMG.  Beginning  with  the  full  order  angular  motion  equations  equipped  with  a 
momentum  exchange  device  such  as  a  MSGCMG  given  by  Wie  (1998,  p.  437) 

sH  +  (dx5H  =  T  (1) 

where  sU(t)eR\  V  s H  =  sHl ,  SH2 ,  7/3  J  the  total  system’ s  angular  momentum  vector 
is  expressed  in  the  body-fixed  control  frame  with  respect  to  an  inertial  frame, 
Text  (/)  e  R3,  VTa,,  =  \Thxt  x ,  Text  2,Thxt  3  J  is  the  external  torque  vector  acting  on  the 

body-fixed  frame  and  to  ( / )  e  M3,  Vto  =  [<Uj,<u2,&»3]r  is  the  angular  velocity  vector  of  the 
body-fixed  frame  with  respect  to  an  inertial  frame. 

The  total  system’s  angular  momentum  vector  includes  both  the  rigid  body’s 
angular  momentum  and  the  MSGCMG’s  angular  momentum.  Therefore 

sH  =  Jrto  +  h  (2) 

where  JT  e  E 3-1:3  is  the  total  moment  of  inertia  of  the  spacecraft  with  the  SGCMG  such 
that 


JT  JB  +  JCmg 


(3) 


where  JB  e  M3*3 


■  J P  + 


1  0 
0  1 
0  0 


0 

0 

mCMGrCMG 


(4) 


represents  the  sum  of  the  second  moment  of  inertia  dyadic  of  the  SGCMG  platform 
JP  e  E3a'3  and  the  second  moment  of  inertia  of  a  point  mass  concentrated  in  the  center  of 
mass  mCMG  of  the  SGCMG  with  respect  to  the  center  of  mass  of  the  spacecraft  rCMG 
assuming  the  SCCMG  is  aligned  with  the  spacecraft’s  third  control  axis.  JCMG  e  E3rf  is 
the  second  moment  of  inertia  dyadic  of  the  MSGCMG  cluster  represented  in  the  body- 
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fixed  frame.  This  represents  the  summation  of  the  second  moment  of  inertia  dyadic  of 
the  MSGCMG  body  and  includes  both  the  contribution  of  its  spinning  inertia  disc  and 
associated  hardware  (Bevilacqua,  Izzo,  and  Valente  2003).  Finally, 

h(/)el3,Vh  =  [/!p/!2,/!,]r  represents  the  total  momentum  vector  of  the  SGCMG 
expressed  in  the  body-fixed  frame.  Differentiation  of  (2)  yields 


s  H  =  JTd)  +  jT(t)  +  h 
=  Jj,6)  +  co  JT(o  +  h 


(5) 


where  of  e  IR3*3  represents  the  skew  symmetric  matrix  with  respect  to  co 


co 


X 


0 

-ru3 

0)2 

c o , 

0 

-cox 

-of 

cox 

0 

(6) 


By  combining  (1),  (2),  and  (5),  we  get 

/I.cb  +  ryxJr7.co  +  h  +  cox(/7.co  +  h)  =  T 'EXT  (V) 

Furthermore,  by  introducing  the  control  torque  vector  generated  on  the  body-fixed 
frame  by  the  MSGCMG  denoted  as  T CMG(t)  e  R3xl ,  VTov/c  =  \_TCMG  VTCMG2,TCMG^  and  a 

control  vector  T(7)  e  U3xl,  VT  =  \Tx,T2, T3]t  to  represent  the  torque  required  by  a  given 
control  law,  (7)  becomes 

J T  co  +  cxf  J  T  co  =  T CMG  +  f f.-xT  (8) 

where 

TCMG  =-(<uxJrco  +  h  +  coxh)  =  T  (9) 

(9)  can  be  solved  for  the  time  derivative  of  the  angular  momentum  vector  as 

h  = -T  -  cox  Jrco  -  co  x  h  (10) 
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g 

▼c 

Figure  4.  Orthonormal  Vectors  of  MSGCMG  Unit  (After  Kurokawa  1998) 

Figure  4  depicts  the  internal  vectors  of  the  MSGCMG  unit.  Three  mutually 
orthogonal  unit  vectors  exist  where  g  is  the  gimbal  vector,  h  is  the  angular  momentum 
vector,  and  c  is  the  torque  vector  where 

3h  ,  .... 

c  =  —  =  gxh  (11) 

do 


The  gimbal  vectors  are  constant  while  the  others  are  dependent  upon  the 
MSGCMG’ s  gimbal  angle  8  .  The  total  angular  momentum  becomes  a  general  function 
of  the  CMG  gimbal  angles  5  and  the  constant  angular  momentum  of  the  unit’s  rotor 
wheel  denoted  by  hw  such  that 


h  =  hw  [0, cosc>, sin  c>]r  (12) 

The  total  output  torque  of  the  MSGCMG  in  the  absence  of  coupling  terms  due  to 
spacecraft  motion  then  can  be  found  by  taking  the  time  derivative  of  (12)  as 

Tcmg  =  -h  =  [0,  K  sin  88,  -hw  cos  88  J  (13) 

By  taking  the  third  body-fixed  axis  to  be  the  yaw  axis  as  depicted  in  Figure  4,  and 
considering  only  the  components  about  this  axis  due  to  compensation  of  the  other  two 
axes  by  the  reactions  of  the  floor,  (8)  and  (9)  simplify  to 

J3dx  =T3-  -hw  cos  88  (14) 
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To  determine  the  steering  logic  for  the  instantaneous  angular  rate  to  command  the 
MSGCMG’s  gimbal  motor,  it  is  thus  necessary  to  solve  (14)  for  the  gimbal  rate  8  given 
the  desired  torque  as  generated  by  a  given  attitude  control  law.  Once  the  gimbal  angular 
position  reaches  ±/r/2 ,  the  MSGCMG  cannot  generate  any  torque  about  the  third  body- 
fixed  axis  and  it  is  therefore  in  a  saturated  singular  configuration. 

2.  Reduced  Order  Dynamics  with  Dual  Rotating  Thrusters  and 
MSGCMG 

The  translation  and  attitude  motion  of  the  simulator  are  governed  by  the  equations 

*BI  =  Vai 

V  =  m~x  'R 

\  Bi  ,n  1vb  A  (15) 

0,  =  o), 

6),  =  J3 1 7) 

where  ,;F(/|  e M2,  VWF  =  \^BFV  BF2 J  represents  the  thrusters’  inputs  limited  to  the  region 
±7t/ 2  with  respect  to  each  face  normal  and  f3(f)el  is  the  attitude  input.  1 RB  , /;  F  and 
T,  are  given  by 


C03 

-S03 

s6*3 

c63 

(16) 


B¥r  =  B¥r  +  B¥T=\_FaC  ^  +FbCab, -Fa  s  +Fbsabf  (17) 

T3=[TcMG+L{-Fasaa-Fbsocb)-\  (18) 

where  s*  =  sin(*),  c»  =  cos(*) . 

The  internal  dynamics  of  the  vectorable  thrusters  are  assumed  to  be  linear 
according  to  the  following  equations 

<=PaJa=  J~Xtab  =  fib,  A  =  J'X  (19) 
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where  J a  and  J h  represent  the  moments  of  inertia  about  each  thruster  rotational  axis 
respectively  and  Ta  e  M ,  Tb  e  R  represent  the  corresponding  thruster  rotation  control 
input. 

The  system’s  state  equation  given  by  (15)  can  be  rewritten  in  control-affine 
system  form  as  (LaValle  2006) 

Nu 

x  =  f  (x)  +  2>,.g,  (x)  =  f  (x)  +  G(x)u,  x  e  (20) 

;=i 

where  Nu  is  the  number  of  controls.  With  M  Vr  representing  a  smooth  Nx  -dimensional 
manifold  defined  be  the  size  of  the  state-vector  and  the  control  vector  to  be  in  Uw" , 
Defining  the  state  vector  xel10  as  \T  =  [x15x2,...,x10]  = 

l>B/,i  •>  rR,,2  ^3’aa->  ab  >  ,  a3  >  Pa  ’  fit,  ] ancl  the  control  vector  ueU5  as 

u'  =[u1,u2,...,u5]  =  [BFa,BFb,TCMG,Ta,Tb],  the  system’s  state  equation  becomes 


X  f  (X)  +  G(x)u  |^Xg  ,  Xy  ,  Xg  ,  Xg  ,  Xjq 


,0j  + 


(x). 


u 


(21) 


where  the  matrix  G,  ( x )  e  R5x5  is  obtained  from  (15)  as 


G,(x)  = 


-m1  [c  x3  c  x4  -  s  x3  s  x4  ] 
-m  1  [c  x3  s  x4  +  s  x3  c  x4  ] 
1  s  x4L 
0 
0 


/?r‘[cx3cx5-sx3sx5]  0 

/?r'[cx3sx5+sx3cx5]  0 

-J, 1  s  x5  L  J, 1 

0  0 

0  0 


0 

0 

0 

j:1 


0 

0 

0 

0 


o  j : 


(22) 


With  the  system  in  the  fonn  of  (20)  given  the  vector  fields  in  (21)  and  (22),  and 
given  that  f(x)  e  M1"  (the  drift  tenn)  and  G(x)el10,s  (the  control  matrix  of  control 
vector  fields)  are  smooth  functions,  it  is  important  to  note  that  it  is  not  necessarily 
possible  to  obtain  zero  velocity  due  to  the  influence  of  the  drift  term.  This  fact  places  the 
system  in  the  unique  subset  of  control-affine  systems  with  drift  and,  as  seen  later,  will 
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call  for  an  additional  requirement  for  detennining  the  controllability  of  the  system. 
Furthermore,  when  studying  controllability  of  systems,  the  literature  to  date  restricts  the 
consideration  to  cases  where  the  control  is  proper.  Having  a  proper  control  implies  that 
the  affine  hull  of  the  control  space  is  equal  to  U^"  or  that  the  smallest  subspace  of  U  is 
equal  to  the  number  of  control  vectors  and  that  it  is  closed  (Sussman  1987;  Sussman 
1990;  Bullo  and  Lewis  2005;  LaValle  2006).  With  a  system  such  as  a  spacecraft  in 
general  or  the  simplified  model  of  the  3 -DoF  simulator  in  particular,  the  use  of  on-off 
cold-gas  thrusters  restrict  the  control  space  to  only  positive  space  with  respect  to  both 
thrust  vectors  leading  to  an  unclosed  set  and  thus  improper  control  space.  In  order  to 
overcome  this  issue,  a  method  which  leverages  the  symmetry  of  the  system  is  used  by 
which  the  controllability  of  the  system  is  studied  by  considering  only  one  virtual  rotating 
thruster  that  is  positioned  a  distance  L  from  the  center  of  mass  with  the  vectored  thrust 
resolved  into  a  y  and  x-component.  In  considering  this  system  perspective,  the  thruster 
combination  now  spans  M2  and  therefore  is  proper  and  is  analogous  to  the  planar  body 
with  variable-direction  force  vector  considered  in  (Lewis  and  Murray  1997;  Bullo  and 
Lewis  2005).  Furthermore,  under  the  assumption  that  the  control  bandwidth  of  the 
thrusters’  rotation  is  much  larger  than  the  control  bandwidth  of  the  system  dynamics,  the 
internal  dynamics  of  the  vectorable  thrusters  can  be  decoupled  from  the  state  and  control 
vectors  for  the  system  yielding  a  thrust  vector  dependent  on  simply  a  commanded  angle. 
Thus  the  system’s  state  vector,  assuming  that  both  thrusters  and  a  momentum  exchange 
device  are  available,  becomes  xr  =  [xj,x2,...,x6]  =  [rBI  vrBI  2,03,vSJ  2,&>3]  e  IR6  and 

the  control  vector  is  nT  =[ul,u2,u3]=[BFl,BF2,Ti]  =  Ui  so  that  the  system’s  state 
equation  becomes 


x  =  f  (x)  +  G(x)u  =  [x4,x5,x6,0,0,0]r  + 


'3x3 


G.W. 


u 


(23) 


where  f  (x)  e  M6  and  G(x)  e  M6*3  are  again  smooth  functions.  The  matrix  Gl  (x)  can  be 
obtained  by  considering  the  relation  of  the  desired  control  vector  to  the  body  centered 
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reference  system,  in  the  two  cases  of  positive  force  needed  in  bx  ( BUl  >  0  )  and  negative 
force  needed  in  bx  (BUX  <0).  In  this  manner,  the  variables  in  (22)  and  (23)  can  be 
defined  as 


ju T  =[BFl,BF2,T3]  =  [~Fa  c  x4 , -Fa  s x4 , Tcmg  ] 
\d  =  -L,Fh=0 

ju7  =  [BFl,BF2,T3]  =  [Fhcx5,Fhsx5,TCMG\ 

\d  =  L,Fa=  0 


(24) 


yielding  the  matrix  in  Gx  (x)  e  M:U3  through  substitution  into  (22)  as 


G,(x) 


m 

cx3 

—m 

sx3 

1  sx3 

m~ 1 

cx3 

o  -dj;1  j~l 


(25) 


When  the  desired  control  input  to  the  system  along  bx  is  zero,  both  thrusters  can 

be  used  to  provide  a  control  force  along  b2 ,  while  a  momentum  exchange  device  provides 
any  required  torque.  In  this  case,  the  control  vector  u  becomes 
uy  =  [uvu2]  =[bF2,T3\  =  U2  and  control  matrix  of  control  fields  becomes  G(x)  e  M6x2in 
(23)  such  that  the  variables  in  (22)  and  (23)  can  be  defined  as 


u T  =[BF2,T}]  =  [Fsa,TCMG\ 

F  =  Fa  =  Fb,a  =  x4  =  -x5  =  ~sign(BU2) 


(26) 


yielding  the  matrix  Gx  (x)  e  K3a"  through  substitution  into  (22)  as 


G,(x) 


-2m  3sx3  0 
2m_Icx3  0 
0  Jf1 


(27) 
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As  will  be  demonstrated  later,  the  momentum  exchange  device  is  not  necessary  to 
ensure  small  time  controllability  for  this  system.  In  considering  this  situation,  which  also 
occurs  when  a  control  moment  gyroscope  is  present  but  is  near  the  singular  conditions 
and  therefore  requires  desaturation,  the  thruster  not  being  used  for  translation  control  can 
be  slewed  to  ±7t/2  depending  on  the  required  torque  compensation  and  fired  to  affect  the 
desired  angular  rate  change.  The  desired  control  input  to  the  system  with  respect  to 
bl  ( BUl  j  can  again  be  used  to  define  the  desired  variables  such  that 


BU,  <  0  -» 
BUX  >  0  -» 


(u7  =[BFl,BF2,T3]  =  [BFl,BF2, 
{ d  =  -L,  x5  =  ±7r/ 2 
|u T  =[BFl,BF2,T}]  =  [BFl,BF2, 
=  L,  x4  =  ±7r/ 2 


^2/to(3  ]  [-Facx4,-Fasx4,Fbdsx5\ 
Fllrofi  ]  =  [Fhcx5,Fhsx5,-Fadsx4] 


(28) 


yielding  the  matrix  G]  (x)  e  K3a!  through  substitution  into  (22)  as 


G,(x) 


III 

cx3 

-m 

sx3 

m 

1  sx3 

m  1 

cx3 

0  -djf 


-[md ) 
( md ) 


sx. 


ex. 


j; 


(29) 


In  case  of  zero  force  requested  along  x  with  only  thrusters  acting,  the  system 
cannot  in  general  provide  the  requested  torque  value. 

A  key  design  consideration  with  this  type  of  control  actuator  configuration  is  that 
with  only  the  use  of  an  on/off  rotating  thruster  to  provide  the  necessary  torque 
compensation,  fine  pointing  can  be  difficult  and  more  fuel  is  required  to  affect  a  desired 
maneuver  involving  both  translation  and  rotation. 
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C.  SMALL-TIME  LOCAL  CONTROLLABILITY 

Before  studying  the  controllability  for  a  nonlinear  control-affine  system  of  the 
form  in  (20),  it  is  important  to  review  several  definitions.  First,  the  set  of  states  reachable 
in  time  at  most  T  is  given  by  /L  (x0,<  T)  by  solutions  of  the  nonlinear  control-affine 
system. 

Definition  1  (Accessibility) 

A  system  is  accessible  from  x0  (the  initial  state)  if  there  exists  r  >  Osuch  that  the 
interior  of  fC  (x0,<  t)  is  not  an  empty  set  for  t  e  ]0,  r]  (Bullo  and  Lewis  2005). 

Definition  2  (Proper  Small-time  Local  Controllability) 

A  system  is  small-time  locally  controllable  (STLC)  from  x0  if  there  exists  r  >  0 
such  that  x0  lies  in  the  interior  of  (x0,<  t )  for  each  t  e  ]0,  r]for  every  proper  control 
set  U  (Bullo  &  Lewis,  2005).  Assuming  that  at  x(0)  =  0  this  can  also  be  seen  under  time 
reversal  as  the  equilibrium  for  the  system  x0  can  be  reached  from  a  neighborhood  in 
small  time  (Sussman  1987;  Sussman  1990). 

Definition  3  (Proper  Control  Set) 

A  control  set  u1  =[uv...,uk]  is  tenned  to  be  proper  if  the  set  satisfies  a  constraint 

ueAi  where  affinely  spans  LJ k  (Sussman  1990;  Bullo  and  Lewis  2005;  LaValle 
2006). 

Definition  4  (Lie  derivative) 

The  Lie  derivative  of  a  smooth  scalar  function  g(x)  e  M  with  respect  to  a  smooth 
vector  field  f  (x)  e  R'v'  is  a  scalar  function  defined  as  (Slotine  1991,  p.  229) 


L{g  =  Vg  f 


r 

1 

/i(x) 

dg_ 

dg 

dx. 

dxN 

1 

ly  x 

_/vr(x)_ 

(30) 
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Definition  5  (Lie  Bracket) 

The  Lie  bracket  of  two  vector  fields  f  (x)  e  R‘V)r  and  g(x)  e  is  a  third  vector 
field  [f.g]  e  M,v'  defined  by  [f,g]  =  Vg  f-Vf  g  ,  where  the  i-th  component  can  be 
expressed  as  (Slotine  1991,  p.  230) 


[f.g],=Z 


'  dXj  '  dx  :  J 


j= 1 


(31) 


Using  Lie  bracketing  methods  that  produce  motions  in  directions  that  do  not  seem 
to  be  allowed  by  the  system  distribution,  sufficient  conditions  can  be  met  to  detennine  a 
system’s  STLC  even  in  the  presence  of  a  drift  vector  as  in  the  equations  of  motion 
developed  above.  These  sufficient  conditions  involve  the  Lie  Algebra  Rank  Condition 
(LARC). 

Definition  6  (Associated  Distribution  A(x) ) 

Given  a  system  as  in  (6),  the  associated  distribution  A(x)  is  defined  as  the  vector 
space  (subspace  of  M,v< )  spanned  by  the  system  vector  fields  figp.-.g^  . 

Definition  7 

The  Lie  algebra  of  the  associated  distribution  ./’(A)  is  defined  to  be  the 

distribution  of  all  independent  vector  fields  that  can  be  obtained  by  applying  subsequent 
Lie  bracket  operations  to  the  system  vector  fields.  Of  note,  no  more  than  Nx  vector 

fields  can  be  produced  (LaValle  2006).  With  dim  [2f  (A)]  <  Nx ,  the  computation  of  the 

elements  of  _Z’  (  A )  ends  either  when  Nx  independent  vector  fields  are  obtained  or  when 
all  subsequent  Lie  brackets  are  vector  fields  of  zeros. 

Definition  8  (Lie  Algebra  Rank  Condition  [LARC]) 

The  Lie  Algebra  Rank  Condition  is  satisfied  at  a  state  x  if  the  rank  of  the  matrix 
obtained  by  concatenating  the  vector  fields  of  the  Lie  algebra  distribution  at  x  is  equal  to 
Nx  (the  number  of  state). 
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For  a  driftless  control-affine  system,  following  the  Chow-Rashevskii  Theorem, 
the  system  is  STLC  if  the  LARC  is  satisfied  (Lewis  and  Murray  1997;  Bulk)  and  Lewis 
2005;  LaValle  2006).  However,  given  a  system  with  drift,  in  order  to  determine  the 
STLC,  the  satisfaction  of  the  LARC  it  is  not  sufficient:  In  addition  to  the  LARC,  it  is 
necessary  to  examine  the  combinations  of  the  vectors  used  to  compose  the  Lie  brackets  of 
the  Lie  algebra.  From  Sussman’s  General  Theorem  on  Controllability,  if  the  LARC  is 
satisfied  and  if  there  are  no  ill  formed  brackets  in^C(A),  then  the  system  is  STLC  from 

its  equilibrium  point  (Sussman  1987).  The  Sussman’s  theorem,  formally  stated  is 
reported  here  below. 

Theorem  1  (Sussman’s  General  Theorem  on  Controllability) 

Consider  a  system  given  by  (20)  and  an  equilibrium  point  p  e  M'Vr'1  such  that 
f  (p)  =  0  .  Assume  £  (A)  satisfies  the  LARC  at  p  .  Furthermore,  assume  that  whenever 
a  potential  Lie  bracket  consists  of  the  drift  vector  f  (x)  appearing  an  odd  number  of 
times  while  gt  (x),...^^  (x)  all  appear  an  even  number  of  times  to  include  zero  times 

(indicating  an  ill  fonned  Lie  bracket),  there  are  sufficient  successive  Lie  brackets  to 
overcome  this  ill  formed  Lie  bracket  to  maintain  LARC.  Then  the  system  is  STLC  from 
p  (Sussman  1987;  Sussman  1990). 


As  it  is  common  in  literature,  an  ill  formed  bracket  is  dubbed  a  “bad”  bracket 
(Sussman  1987;  Sussman  1990;  Lewis  and  Murray,  1997,  Bulk)  and  Lewis,  2005; 
LaValle  2006).  Conversely,  if  a  bracket  is  not  “bad,”  it  is  termed  “good”.  As  an 
example,  for  a  system  with  a  drift  vector  and  two  control  vectors,  the  bracket 
[f ,  [g1?  gj  ]]  is  bad,  as  the  drift  vector  occurs  only  once  while  the  first  control  vector 


appears  twice  and  the  second  control  vector  appears  zero  times.  Similarly,  the  bracket 


is  good  as  the  first  control  vector  appears  only  once.  Therefore,  it  can  be 


summarized  that  if  the  rank  of  the  Lie  algebra  of  a  control-affine  system  with  drift  is 
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equal  to  the  number  of  states  and  there  exist  sufficient  “good”  brackets  to  overcome  the 
“bad”  brackets  to  reach  the  required  LARC  rank,  then  the  system  is  small  time  locally 
controllable. 

D.  SMALL-TIME  LOCAL  CONTROLLABILITY  CONSIDERATIONS  FOR 

THE  3-DOF  SPACECRAFT  SIMULATOR 

The  concept  of  small  time  local  controllability  is  better  suitable  than  the  one  of 
accessibility  for  the  problem  of  spacecraft  rendezvous  and  docking,  as  a  spacecraft  is 
required  to  move  in  any  directions  in  a  small  interval  of  time  dependent  on  the  control 
actuator  capabilities  (e.g.,  to  avoid  obstacles).  The  finite  time  r  can  be  arbitrary  if  the 
control  input  is  taken  to  be  unbounded  and  proper  (Sussman  1990;  Bullo  and  Lewis  2005; 
LaValle  2006). 

While  no  theory  yet  exists  for  the  study  of  the  general  controllability  for  a  non¬ 
linear  system,  the  STLC  from  an  equilibrium  condition  can  be  studied  by  employing 
Sussman’s  theorem.  For  the  case  of  spacecraft  motion,  in  order  to  apply  Sussman’s 
theorem,  we  hypothesize  that  the  spacecraft  is  moving  from  an  initial  condition  with 
velocity  close  to  zero  (relative  to  the  origin  of  an  orbiting  reference  frame). 

In  applying  Sussman’s  General  Theorem  on  Controllability  to  the  reduced  system 
equations  of  motion  presented  in  (23)  with  G,  (x)  given  in  (25),  the  Lie  algebra  evaluates 
to 


(A)  =  span  {gj ,  g2,  g3  [f ,  gi  ] ,  [f ,  g2  ] ,  [f ,  g3  ]}  (32) 

so  that  dim[zf (A)]  =  Nx  =  6.  In  order  to  verify  that  this  is  the  minimum  number  of 

actuators  required  to  ensure  STLC,  the  Lie  algebra  is  reinvestigated  for  each  possible 
combination  of  controls.  Appendix  A  includes  a  developed  MATLAB®  function  for 
determining  the  Lie  algebra  for  an  arbitrary  nonlinear  system.  The  resulting  analysis,  as 
summarized  in  Table  2,  demonstrates  that  the  system  is  STLC  from  the  systems 
equilibrium  point  at  f  (p)  =  0  given  either  two  rotating  thrusters  in  complementary  semi¬ 
circle  planes  or  fixed  thrusters  on  opposing  faces  providing  a  nonnal  force  vector  to  the 
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face  in  opposing  directions  and  a  momentum  exchange  device  about  the  center  of  mass. 
For  instance,  in  considering  the  case  of  control  inputs  ur  =  [i{l,u2\  =  [BFy,T:]eU2  ,(23) 
given  (27)  becomes 

x  =  f(x)  +  g,(x)n1+g2(x)w2 

T  T  T  (33) 

=  [x4,x5,x6, 0,0,0]  +^0,0,0,2m~1sx3,2m~1cx3,0]  ux  +[0,0,0,0,0, J3l J  u2 

with  the  equilibrium  point  p  such  that  f(p)  =  0  is  p  =  [x15x2,x3,0,0,0]7  .  The  ./’(A)  is 
formed  by  considering  the  associated  distribution  A(x)  and  successive  Lie  brackets  as 


f, 

gl> 

g2 

[f,8i]’ 

[gpg2]. 

[f’§2] 

[f»[gpg2]]» 

[f’[f’g2]]’ 

[g!,[f,gl]]. 

[gi>[gi>g2]]» 

[gl,[f5g2]]? 

[g2>[f’gl]]’ 

[  g  2 5  [  g  1  **  g  2  ]  ]  ^ 

[g2’[f’g2]]’ 

[f,[f,[f,g2]]]. 

[f,[f,[gl,g2]]]. 

[f,[f’[f’g2]]]’ 

The  sequence  can  first  be  reduced  by  considering  any  “bad”  brackets  in  which  the 
drift  vector  appears  an  odd  number  of  times  and  the  control  vector  fields  each  appear  an 
even  number  of  times  to  include  zero.  In  this  manner  the  Lie  brackets  [gt , [f - g3  ]] 

and  [g2 ,  [f ,  g2  ]]  can  be  disregarded. 

By  evaluating  each  remaining  Lie  bracket  at  the  equilibrium  point  p ,  the  linearly 
independent  vector  fields  can  be  found  as 
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gi  =  [03jci ,  — m  'sx3,m  'cXjj-Jj'zJ 
g2=[05xl^3"1T 

[f ,gi ]  =  Vgi -f-Vf.g,  ^\_m~lsx3,-m~lcx3,04xlJ 

[f ,  g2  ]  =  Vg2  •  f  -  Vf  •  g2  =  [oM ,  -  j;1 , 03xl  ]r 

[gi» [f » g2 ]]  =  v  [f . g2  ]  •  gi  -  Vgj-  [f , g2 ]  =  [o3li, - (j mJ3 )~‘  cx3 - (mJ3 )“‘  sx3, 0 

[f ,  [gi ,  [f ,  g2  ]]]  =  V  [gl ,  [f ,  g2  ]]  •  f  -  Vf  •  [g, ,  [f ,  g2  ]]  =  [2L  ( mJ3  )~l  cx3 , 2 L  ( mJ3 )“'  5X3 , 04x, 

(34) 

Therefore,  the  Lie  algebra  comprised  of  these  vector  fields  is 


(a)  =  spcin  { gi ,  g2  9  [f ,  gj  ] ,  [f ,  g2  ] ,  [g3 ,  [f ,  g2  ]] ,  [f ,  [gi ,  [f ,  g2  ]]]}  (35) 


yielding  dim[Z’(A)]  =  Nx  =  6,  and  therefore  the  system  is  small  time  locally 
controllable. 


Control  Inputs 

Thruster  Positions 

dim(£(A)) 

Controllability 

bF, 

o 

ll 

if 

II 

a 

8 

2 

Inaccessible 

bf2 

«a  =~ab  =  W2 

5 

Inaccessible 

T 

1 CMG 

NA 

2 

Inaccessible 

Bp  T 

1  1  ’ ±  CMG 

aa=ah=0 

6 

STLC 

(BF  T  )or(BF  bF  ) 

\  1  2  ’ 1  CMG  )Ur  \  1  2  9  1  2/ rot  3  ) 

aa=~ah=+x/2 

6 

STLC 

Bf  BJ7 

1  P  1  2 

\a\<nll\ab\<n/2 

6 

STLC 

Table  2.  STLC  Analysis  for  the  3-DoF  Spacecraft  Simulator 

E.  NAVIGATION  AND  CONTROL  OF  THE  3-DOF  SPACECRAFT 
SIMULATOR 


In  the  current  research,  the  assumption  is  made  that  the  spacecraft  simulator  is 
maneuvering  in  the  proximity  of  an  attitude  stabilized  target  spacecraft  and  that  this 
spacecraft  follows  a  Keplarian  orbit.  Furthermore,  the  proximity  navigation  maneuvers 
are  considered  to  be  fast  with  respect  to  the  orbital  period.  A  pseudo-GPS  inertial 
measurement  system  by  Metris,  Inc.  (iGPS)  is  used  to  fix  the  ICS  in  the  laboratory  setting 
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for  the  development  of  the  state  estimation  algorithm  and  control  commands.  The  X-axis 
is  taken  to  be  the  vector  between  the  two  iGPS  transmitters  with  the  Y  and  Z  axes  forming 
a  right  triad  through  the  origin  of  a  reference  system  located  at  the  closest  corner  of  the 
epoxy  floor  to  the  first  iGPS  transmitter.  Navigation  is  provided  by  fusing  of  the 
magnetometer  data  and  fiber  optic  gyro  through  a  discrete  Kalman  filter  to  provide 
attitude  estimation  and  through  the  use  of  a  linear  quadratic  estimator  to  estimate  the 
translation  velocities  given  inertial  position  measurements.  Control  is  accomplished 
through  the  combination  of  a  state  feedback  linearized  controller,  a  linear  quadratic 
regulator,  Schmitt  trigger  logic  and  Pulse  Width  Modulation  using  the  minimal  control 
actuator  configuration  of  the  3 -DoF  spacecraft  simulator.  Figure  5  reports  a  block 
diagram  representation  of  the  control  system. 


Figure  5.  Block  diagram  of  the  control  system  of  the  3-DoF  spacecraft  simulator, 

(vi,  V2,  V3  control  inputs  to  the  system) 
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1.  Navigation  Using  Inertial  Measurements  with  Kalman  Filter  and 
Linear  Quadratic  Estimator 

In  the  presence  of  the  high  accuracy,  low  noise,  high  bandwidth  iGPS  sensor  with 
position  accuracy  to  within  5.4  mm  with  a  standard  deviation  of  3.6  mm  and 
asynchronous  measurement  availability  with  a  nominal  frequency  of  40  Hz,  a  full-order 
linear  quadratic  estimator  with  respect  to  the  translation  states  is  implemented  to 
demonstrate  the  capability  to  estimate  the  inertial  velocities  in  the  absence  of 
accelerometers.  Additionally,  due  to  the  effect  of  noise  and  drift  rate  in  the  fiber-optic 
gyro,  a  discrete -time  linear  Kalman  filter  is  employed  to  fuse  the  data  from  the 
magnetometer  and  the  gyro.  Both  the  gyro  and  magnetometer  are  capable  of  providing 
new  measurements  asynchronously  at  100  Hz.  Experimental  results  are  presented  in 
subsequent  sections  to  demonstrate  the  filter’s  effectiveness. 


a.  Attitude  Discrete-time  Kalman  Filter 

With  the  attitude  rate  being  directly  measured,  the  measurement  process 
can  be  modeled  in  state-space  equation  form  as: 


(36) 


z  =  n,=[l  0]  V  +%m  (37) 

'  H  'L  8  - 

where  co •  is  the  measured  gyro  rate,  fi  is  the  gyro  drift  rate,  ricog  and  are  the 
associated  gyro  output  measurement  noise  and  the  drift  rate  noise  respectively.  y/m  is  the 
measured  angle  from  the  magnetometer,  and  rji//m  is  the  associated  magnetometer  output 
measurement  noise.  It  is  assumed  that  n,  rtRa  and  are  zero-mean  Gaussian  white- 
noise  processes  with  variances  given  by  o)0g,  cl,,  and  a\,m  respectively.  Introducing 
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the  state  variables  xr  =[V,/?g],  control  variables  u=cog,  and  error  variables 
w7  =  ,  t]/jg ]  and  v  =  rji/m ,  (36)  and  (37)  can  be  expressed  compactly  in  matrix  form  as 


x(t)  =  A(t)x(t)  +  B(t)u(t)  +  G(t)w(t) 

(38) 

z(t)  =  Hx{t)  +  v(t) 

(39) 

In  assuming  a  constant  sampling  interval  At  in  the  gyro  output,  the  system 
equation  (38)  and  observation  equations  (39)  can  be  discretized  and  rewritten  as 

x*+i=®*x*+r*»*+T*wt  (40) 

za-  =Hk*k+vk  (41) 


where 


=  e^' 


1  -At 
0  1 


(42) 


and 


(43) 


The  process  noise  covariance  matrix  used  in  the  propagation  of  the 
estimation  error  covariance  given  by  (Gelb  1974;  Crassidis  and  Junkins  2004) 

(44) 


can  be  properly  numerically  estimated  given  a  sufficiently  small  sampling  interval  by 
following  the  numerical  solution  by  van  Loan  (Crassidis  and  Junkins  2004).  First,  the 
following  2 n  x  2 n  matrix  is  fonned: 


-A  GQGr 
0  At 


At 


(45) 
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where  At  is  the  constant  sampling  interval,  A  and  G  are  the  constant  continuous-time 
state  matrix  and  error  distribution  matrix  given  in  (38),  and  Q  is  the  constant  continuous¬ 
time  process  noise  covariance  matrix 


£>  =  £{w(0wr(0} 


The  matrix  exponential  of  (45)  is  then  computed  by 


X 

%  1 

0 

% 

0 

_ 

(46) 


(47) 


where  ®A.  is  the  state  transition  matrix  from  (42)  and  =  ( YkQk  ft  f  ■  Therefore,  the 
discrete-time  process  noise  covariance  is 


<\={vkQJk)T 


l/3^gAt3  +  (Tg2At  -l/2rr2gAr 
-l/2o>gAr  a\gAt 


(48) 


The  discrete -time  measurement  noise  covariance  is 

rk=E{vyk}=(jim  (49) 

Given  the  fdter  model  as  expressed  in  (36)  and  (37),  the  estimated  states 
and  error  covariance  are  initialized  where  this  initial  error  covariance  is  given 
by  PQ  =  A’  |  \(/0)\r(/0)| .  If  a  measurement  is  given  at  the  initial  time,  then  the  state  and 
covariance  are  updated  using  the  Kalman  gain  formula 

(50) 


where  P'k  is  the  a  priori  error  covariance  matrix  and  is  equal  to  P0 .  The  updated  or  a 
posteriori  estimates  are  determined  by 
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(51) 


rk  =  xk  +  Kk  [zk  -Hkxk  J 
p;=[hX2  -KkHk]p- 


where  again  with  a  measurement  given  at  the  initial  time,  the  a  priori  state  xk  is  equal  to 
x0 .  The  state  estimate  and  covariance  are  propagated  to  the  next  time  step  using 


xi+i  =  ®A+rVb 


(52) 


If  a  measurement  is  not  given  at  the  initial  time  step,  or  any  time  step 
during  the  process,  the  estimate  and  covariance  are  propagated  to  the  next  available 
measurement  point  using  (52). 


b.  Translation  Linear  Quadratic  Estimator 

With  the  measured  translation  state  from  the  iGPS  sensor  being  given  by 


z  = 


10  0  0 
0100 

c 


(53) 


the  dynamics  of  a  full-order  state  estimator  is  described  by  the  equation 

x  =  x-  x  =  [^lx  +  J8u]-[ylx  +  f?u  +  Llqe  (z - Cx) J 
=  A(x-x)-  Llqe(Cx-Cx)  (54) 

—  ^A-LlqeC^x 


where 


Ax  +  Bu  :  linearized  plant  dynamics 

Ax  +  Bu  :  system  model 

Lloe  :  linear  quadratic  estimator  gain  matrix 

Cx  :  predicted  measurement 
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The  observer  gain  matrix  Lloe  can  be  solved  using  standard  linear 
quadratic  estimator  methods  as  (Bryson  1993) 

Llqe  =  PCT  Rtx  (55) 

where  P  is  the  solution  to  the  algebraic  Riccati  equation 

AP  +  PAt  -  PCTRT  lCP  +  Qt  =  0  (56) 


and  Qt  and  Rr  are  the  associated  weighting  matrices  with  respect  to  the  translational 
degree  of  freedom  defined  as 


Qt  =  diag  (l  / 1 Armax  f ,  1  / 1| Armax  f ,  1  / 1| Avmax  f  ,  1  / 1| Avmax  f ) 
RT=diag(\/(FmJ,\/(FmJ) 


(57) 


where  Ar  I ,  AvI  are  taken  to  be  the  maximum  allowed  errors  between  the  current 
and  estimated  translational  states  and  Fmax  is  the  maximum  possible  imparted  force  from 
the  thrusters. 

Table  3  lists  the  values  of  the  attitude  Kalman  filter  and  translation  state 
observer  used  for  the  experimental  tests. 
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2.  Smooth  Feedback  Control  via  State  Feedback  Linearization  and 
Linear  Quadratic  Regulation 

Considering  a  Multi-Input  Multi-Output  (MIMO)  nonlinear  system  in  control- 
affine  form,  the  state  feedback  linearization  problem  of  nonlinear  systems  can  be  stated 
as  follows:  obtain  a  proper  state  transformation 

z  =  ®(x)  where  z  e  R ,Vt  (58) 


and  a  static  feedback  control  law 


u  =  a(x)  +  /?(x)v  where  velR^"  (59) 

such  that  the  closed-loop  system  in  the  new  coordinates  and  controls  become 

v  (60) 

x=0-1(z) 

is  both  linear  and  controllable.  The  necessary  conditions  for  a  MIMO  system  to  be 
considered  for  input-output  linearization  are  that  the  system  must  be  square  or  Nu  =  Ny 
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where  Nu  is  defined  as  above  to  be  the  number  of  control  inputs  and  Nv  is  the  number 
of  outputs  for  a  system  of  the  expanded  form  (Isidori  1989;  Slotine  1990) 


x  =  f  (x)  +  G(x)u 

Ny 

y  =  ZA(x)  =  h(x) 

i= 1 


(61) 


The  input-output  linearization  is  determined  by  differentiating  the  outputs  y.  in 
(61)  until  the  inputs  appear.  Following  the  method  outlined  in  Slotine  (1990,  chap.  6)  by 
which  the  assumption  is  made  that  the  partial  relative  degree  r.  is  the  smallest  integer 

such  that  at  least  one  of  the  inputs  appears  in  y.n> ,  then 


*(")=V'A,+£i„v,",A-W“/  <62> 

7=1 

with  the  restriction  that  Lg  Z,f';_1/zf  (x)  ^  0  for  at  least  one  j  in  a  neighborhood  of  the 
equilibrium  point  x0 .  Letting 


E(x) 


h,Lr\(*)  ■■ 

•  V 

LtlL,r’-%(x)  ■■ 

A,,  1 

AA  A  (x) 

•  L  L. 

%N„  1 

rN ..  “I  . 


(63) 


so  that  (63)  is  in  the  form 


1 

_ 1 

J5- 

(a*2  ) 

y2 

— 

4WM 

i 

_ i 

1 

s- 

_ 1 

(64) 


the  decoupling  control  law  can  be  found  where  the  NyxNv  matrix  L(x)  is  invertible 

over  the  finite  neighborhood  of  the  equilibrium  point  for  the  system  as 
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(65) 


v,- V'^iO) 

v2  -L/2h2  (x) 

VNy  ~L{  hNy  (x) 

With  the  above  stated  equations  for  the  simulator  dynamics  in  (23)  given  Gl  (x) 
as  defined  in  (25),  if  we  choose 

h(x)  =  [rs/^f  (66) 

the  state  transformation  can  be  chosen  as 

zT  =  [hl(x),h2(x),h3(x),Lfh1(x),Lth2(x),Lfh3(x)]  =  [rBI,02,\BI,o)3]  (67) 

where  z(t)  e  M6,  Vz  =  [zj,z2,...,z6]r  are  new  state  variables,  and  the  system  in  (23)  is 
transfonned  into 

z  =  z4,z5,z6,/n_1  (cz3  bF{  -sz3  BF2^,m  l  (sz3  BFX  +cz3  BF2^J,J3lT3  (68) 

The  dynamics  given  by  (23)  considering  the  switching  logic  described  in  (24), 
(26)  and  (28)  can  now  be  transformed  using  (68)  and  the  state  feedback  control  law 

[eF,7';]  =  £(x)  1  (v-b)  (69) 

into  a  linear  system 


(70) 


where 


b  =  L)hx  (x),Z'f2/i2  (x),L'f7i3  (x)J 


(71) 
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and  E(x)  given  by  (63)  with  equivalent  inputs  v(/)  e  M3,  Vv  =  [v,,v2,v3]r  and  relative 

degree  of  the  system  at  the  equilibrium  point  x0  is  (r15r2,r3)  =  (2,2,2)  .  Therefore  the 

total  relative  degree  of  the  system  at  the  equilibrium  point,  which  is  defined  as  the  sum  of 
the  relative  degree  of  the  system,  is  six.  Given  that  the  total  relative  degree  of  the  system 
is  equal  to  the  number  of  states,  the  nonlinear  system  can  be  exactly  linearized  by  state 
feedback  and  with  the  equivalent  inputs  v;. ,  both  stabilization  and  tracking  can  be 

achieved  for  the  system  without  concern  for  the  stability  of  the  internal  dynamics  (Slotine 
1990). 

One  of  the  noted  limitations  of  a  feedback  linearized  based  control  system  is  the 
reliance  on  a  fully  measured  state  vector  (Slotine  1990).  This  limitation  can  be  overcome 
through  the  employment  of  proper  state  estimation.  HIL  experimentation  on  SRL’s 
second  generation  robotic  spacecraft  simulator  using  these  navigation  algorithms 
combined  with  the  state  feedback  linearized  controller  as  described  above  coupled  with  a 
linear  quadratic  regulator  to  ensure  the  poles  of  (70)  lie  in  the  open  left  half  plane 
demonstrate  satisfactory  results  as  reported  in  the  following  section. 

a.  Feedback  Linearized  Control  Law  with  CMG  Rotational  Control 
and  Thruster  Translational  Control 

By  applying  (69)  to  the  dynamics  in  (23)  given  Gl  (x)  as  defined  in  (25) 
where  the  system  is  taken  to  be  observable  in  the  state  vector  y  =  [r5/,t?3]r  =  [v, ,x2,x,]7 
and  by  using  thruster  b  for  translational  control  (i.e.,  for  the  case  fiUl  >  0  where 
BUl  =  v,  c  0-,  +  v2  s  03  and  BU2  =  -v,  s03  +  v2c03 ),  the  feedback  linearized  control  law  is 

ur  =  [  BFl ,  bF2  ,  T3  ]  =  [m  BUl ,  m  BU2 ,  ml  BU2  +  J3v3  ]  (72) 

which  is  valid  for  all  x  in  a  neighborhood  of  the  equilibrium  point  x0 .  Similarly,  the 
feedback  linearized  control  law  when  BUl  <  0  (thruster  a  is  providing  translation  control) 

ur  =  [ BF1 ,  bF2 ,T3]  =  [m  BUl , m  BU 2 , -mL  BU 2  +  J3v3 ]  (73) 
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Finally,  when  BU3  =  0  (both  thrusters  used  for  translational  control)  given 
Gj  (x)  as  defined  in  (27)  is 

ur  =  [  bF2 , 7)  ]  =  [m  bU2  / 2 ,  J3v3  ]  (74) 

b.  Feedback  Linearized  Control  Law  for  Thruster  Roto- 
Translational  Control 

As  mentioned  previously,  by  considering  a  momentum  exchange  device 
for  rotational  control,  momentum  storage  must  be  managed.  For  a  control  moment 
gyroscope  based  moment  exchange  device,  desaturation  is  necessary  near  gimbal  angles 
of  rtf  2 .  In  this  region,  due  to  the  mathematical  singularity  that  exists,  very  little  torque 
can  be  exchanged  with  the  vehicle  and  thus  it  is  essentially  ineffective  as  an  actuator.  To 
accommodate  these  regions  of  desaturation,  logic  can  be  easily  employed  to  define 
controller  modes  as  follows:  If  the  MSGCMG  is  being  used  as  a  control  input  and  if  the 
gimbal  angle  of  the  MSGCMG  is  greater  than  75  degrees,  the  controller  mode  is  switched 
from  normal  operation  mode  to  desaturation  mode  and  the  gimbal  angle  rate  is  directly 
commanded  to  bring  the  gimbal  angle  to  a  zero  degree  nominal  position  while  the 
thruster  not  being  directly  used  for  translational  control  is  slewed  as  appropriate  to 
provide  torque  compensation.  In  these  situations,  the  feedback  linearizing  control  law  for 
the  system  dynamics  in  (23)  given  Gj  (x)  as  defined  in  (29)  where  thruster  b  is  providing 

translational  control  ( BUX  >  0  ),  and  thruster  a  is  providing  the  requisite  torque  is 

ur  =  [  BF1 ,  bF2  ,  T3  ]  =  \m  BUX  ,{mL  BU2  -  J3v3 )  flL  ,(mL  BU2  +  J3v3 )  /l\  (75) 

Similarly,  the  feedback  linearizing  control  law  for  the  system  assuming 
thruster  a  is  providing  translational  control  ( BUl  <  0)  while  thruster  b  provides  the 
requisite  torque  is 

ur  =  [  BF1 ,  bF2  ,  T3  ]  =  [m  %  ,{mL  BU2  +  J3v3 )  /lL  ,(mL  BU2  -  J3v3 )  /l\  (76) 
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3.  Determination  of  the  Thruster  Angles,  Forces  and  CMG  Gimbal 
Rates 

In  either  mode  of  operation,  the  pertinent  decoupling  control  laws  are  used  to 
determine  the  commanded  angle  for  the  thrusters  and  whether  or  not  to  open  or  close  the 
solenoid  for  the  thruster.  For  example,  if  BUl  >  0  ,  (72)  or  (75)  can  be  used  to  determine 
the  angle  to  command  thruster  b  as 

ab=\mrl(BF2lBFy)  (77) 

and  the  requisite  thrust  as 

Fh=^BF;  +  BF;  (78) 

If  the  MSGCMG  is  being  used,  the  requisite  torque  commanded  to  the  CMG  is 
taken  directly  from  (72).  In  the  nonnal  operation  mode,  with  the  commanded  angle  for 
thruster  one  not  pertinent,  it  can  be  commanded  to  zero  without  affecting  control  of  the 
system.  Similarly,  if  BU \  <  0 ,  (73)  or  (76)  can  be  used  to  detennine  the  angle  to 

command  thruster  b  and  the  requisite  thrust  analogous  to  (77)  and  (78).  The  requisite 
torque  commanded  to  the  CMG  is  similarly  taken  directly  from  (73).  The  required  CMG 
torques  can  be  used  to  determine  the  gimbal  rate  SCMG  to  command  the  MSGCMG  by 
solving  (14)  given  a  desired  body-fixed  torque  7)  such  that 

& cmg  =  — ^3  /( hw  cos  SCMG )  (79) 

where  hw  is  the  constant  angular  momentum  of  the  rotor  wheel  and  SCMG  is  the  current 
angular  displacement  of  the  wheel’s  rotational  axis  with  respect  to  the  horizontal. 

If  the  momentum  exchange  device  is  no  longer  available  and  BU{  >  0  ,  the  thruster 
angle  commands  and  required  thrust  value  for  the  opposing  thruster  can  be  detennined  by 
using  (75)  as 

aa  --k/2  sign  ( mL  BU 2  +  J3v3  j  (80) 
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and 


Fa  =  -sign(aa )  ( mL  BU  2  +  J3v3 )  /  L  (81) 

given  7)  =  -Fa  s  aaL  .  Likewise,  the  thruster  angle  commands  and  required  thrust  value 
for  the  opposing  thruster  given  BUi  <  0  can  be  detennined  by  using  (76)  as 

ab  =  jt/2  sign  ( mL  BU2  -  J3v3  j  (82) 

and 

Fh  =  sign{ccb)[mL  BU 2  -  J3v3 )/ L  (83) 

given  T3  =  -Fb  sin  abL  . 


4.  Linear  Quadratic  Regulator  Design 

In  order  to  determine  the  linear  feedback  gains  used  to  compute  the  requisite 
equivalent  inputs  v.  to  regulate  the  three  degrees  of  freedom,  so  that 


lim  zi(t)  =  r  (t)  =  r 

f-»oo  ’  \  /  J 

lim  z4(t)  =  vBU  (t)  =  vref. 


\\mz2(t)  =  rB2(t)  =  r 

t — ^oo  ’  x  7  J 

lim  z5{t)  =  vBI2{t)  =  vref, 

t->  00  ’  v  J 


lim  z3{t)  =  63(t)  =  Oref 
^  (84) 

limz6(0  =  0),  (/)  =  ojref 


a  standard  linear  quadratic  regulator  is  employed  where  the  state-feedback  law 
v  =  -Kz  minimizes  the  quadratic  cost  function 


/(v)  =  ^(zr  Qz  +  yT  Rv^dt  (85) 

0 


subject  to  the  feedback  linearized  state-dynamics  of  the  system  given  in  (70)  .  Given  the 
relation  between  the  linearized  state  and  true  state  of  the  system,  the  corresponding  gain 
matrices  R  and  Q  in  (85)  are  chosen  to  minimize  the  appropriate  control  and  state  errors 
as 
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Q  =  diag 


vl/ 1 ^ f  ,  1  / 1 Av^ ||2 , 1  / ( A a>m )2 y 


R  =  diag  (l  /  (Fm  f ,  1  /  (F„  f ,  1  /  (r^  )! ) 


(86) 


where  TCMG  max  is  taken  to  be  the  maximum  possible  imparted  force  and  torques  from  the 
thrusters  and  MSGCMG. 

Given  the  use  of  discrete  cold-gas  thrusters  in  the  system  for  translational  control 
throughout  a  commanded  maneuver  and  rotational  control  when  the  continuously  acting 
momentum  exchange  device  is  unavailable,  Schmitt  trigger  switching  logic  is  imposed. 
Schmitt  triggers  have  the  unique  advantage  of  reducing  undesirable  chattering  and 
subsequent  propellant  waste  nearby  the  reference  state  through  an  output-versus-input 
logic  that  imposes  a  dead  zone  and  hysteresis  to  the  phase  space  as  shown  in  Figure  8. 


Figure  6.  Schmitt  Trigger  Characteristics  with  Design  Parameters  Considering  rB1 , 

Coordinate  Control  Logic 


Three  separate  Schmitt  triggers  are  used  with  the  design  parameters  of  the  Schmitt 
trigger  shown  in  Figure  8  (as  demonstrated  for  the  rBI ,  coordinate  control  logic).  In  the 

case  of  the  two  translational  DoF  Schmitt  triggers,  the  parameters  are  chosen  such  that 


£on  0)  =  Krrdb  +  K,VL 
£off  (  r  )=  Kr  rdb  ~ 


(87) 
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where  vL  =  FmaxAt /2m  .  rdh  is  a  free  parameter  that  is  constrained  by  mission 
requirements.  vout  is  chosen  such  that  the  maximum  control  command  from  the 
decoupling  control  law  yields  a  value  less  than  or  equal  to  Fnax  for  the  translational 
thruster. 

In  the  case  of  the  rotational  DoF  Schmitt  trigger  when  the  momentum  exchange 
device  is  unavailable,  the  parameters  are  chosen  such  that 


8 on  Kg0db  +  KaO)L 

8 off  ~  Kg0db  —  Kwcol 


(88) 


where  coL  =  Fmax  LAt / 2  J3  .  0db  is  a  free  parameter  that  is  again  constrained  by  mission 

requirements.  For  both  modes  of  operation  (i.e.  with  or  without  a  momentum  exchange 
device),  (72)  through  (74)  can  be  used  to  detennine  that 

v.  =  v,  =  F  /42m  (89) 

l,max  2, max  max/  v  \  / 


and  when  the  thrusters  are  used  for  rotational  control 

V3,max  —  FmaJ-'/j^ 


(90) 


When  the  momentum  exchange  device  is  available,  the  desired  torque,  as 
determined  by  the  LQR  control  law  described  above,  is  passed  directly  through  the 
Schmitt  trigger  to  the  decoupling  control  law  to  determine  the  required  gimbal  rate 
command  to  the  MSGCMG. 

The  three  Schmitt  trigger  blocks  output  the  requested  control  inputs  along  the 
I  frame.  The  appropriate  feedback  linearizing  control  law  is  then  used  to  transform 
these  control  inputs  into  requested  thrust,  thruster  angle  and  MSGCMG  gimbal  rate  along 
the  2>  frame.  From  these,  a  vector  of  specific  actuator  commands  are  formed  such  that 

(91) 
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Each  thruster  command  is  normalized  with  respect  to  Fmax  and  then  fed  with  its 

corresponding  commanded  angle  into  separate  Pulse  Width  Modulation  (PWM)  blocks. 
Each  PWM  block  is  then  used  to  obtain  an  approximately  linear  duty  cycle  from  on-off 
actuators  by  modulating  the  opening  time  of  the  solenoid  valves  (Wie  1998,  p.  455). 
Additionally,  due  to  the  linkage  between  the  thruster  command  and  the  thruster  angle,  the 
thruster  firing  sequence  is  held  until  the  actual  thruster  angle  is  within  a  tolerance  of  the 
commanded  thruster  angle.  Furthermore,  in  order  to  reduce  over-controlling  the  system, 
the  LQR,  Schmitt  trigger  logic  and  decoupling  control  algorithm  are  run  at  the  PWM 
bandwidth  of  8.33  Hz.  From  each  PWM,  digital  outputs  (either  zero  or  one)  command 
the  two  thrusters  while  the  corresponding  angle  is  sent  via  RS-232  to  the  appropriate 
thruster  gimbal  motor. 


A6 

max 

1.8  x  10‘2  rad 

A  <x> 

max 

1.8  x  10'2  rad-s'1 

T 

±  CMG,  max 

0.668  Nm 

Kr  =  Klqr(\,\)  =  Klqr(2,2) 

15.9 

Kv  =  Klqr  (1, 4)  =  Klqr  (2, 5) 

84.54  s 

K0  =  Klqr  (3, 3) 

1.39 

Ka  =  Klqr  (3, 6) 

1.75  s 

r*b 

1  O'2  m 

VL 

3.05  x  10'5  m-s'1 

odb 

1.8  x  10'2  rad 

a>L 

1.8  x  10'2  rad-s'1 

£ on  (f) 

1.61  x  10'1  m 

£ofAr) 

1.56  x  10"1  m 

£on(&) 

2.47  x  10"2  rad 

£off(0) 

2.37  x  10"2  rad 

PWM  min  pulse  width 

10'2  s 

PWM  sample  time 

1.2  x  10'1  s 

Table  4.  Values  of  the  Control  Parameters 
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Table  3  and  Table  4  list  the  values  of  the  control  parameters  used  for  the 
experimental  tests  reported  in  the  following  section.  In  particular,  ||Avmax||  is  chosen 

based  typical  maximum  relative  velocities  during  rendezvous  scenarios  while  A #max  is 
taken  to  be  1  degree  and  A  n>max  is  chosen  to  be  1  degrees/sec  which  correspond  to  typical 
slew  rate  requirements  for  small  satellites  (Roser  and  Schedoni  1997;  Lappas,  Steyn,  and 
Underwood  2002).  The  minimum  opening  time  of  the  PWM  was  based  on  experimental 
results  for  the  installed  solenoid  valves  reported  in  (Lugini  and  Romano  2009). 

F.  EXPERIMENTAL  RESULTS 

The  navigation  and  control  algorithms  introduced  above  were  coded  in 
MATLAB®-Simulink®  and  run  in  real  time  using  MATLAB  XPC  Target™  embedded 
on  the  SRL’s  second  generation  spacecraft  simulator’s  on-board  PC-104.  Two 
experimental  tests  are  presented  to  demonstrate  the  effectiveness  of  the  designed  control 
system.  The  scenario  presented  represents  a  potential  real-world  autonomous  proximity 
operation  mission  where  a  small  spacecraft  is  tasked  with  perfonning  a  full  360-degree 
circle  around  another  spacecraft  for  the  purpose  of  inspection  or  pre-docking.  These 
experimental  tests  validate  the  navigation  and  control  approach  and  furthermore 
demonstrate  the  capability  of  the  robotic  spacecraft  simulator  testbed. 

1,  Autonomous  Proximity  Maneuver  Using  Vectorable  Thrusters  and 
MSGCMG  Along  a  Closed  Circular  Path 

Figure  7,  Figure  8,  and  Figure  9  report  the  results  of  an  autonomous  proximity 
maneuver  along  a  closed  circular  trajectory  of  NPS  SRL’s  second  generation  robotic 
spacecraft  simulator  using  its  vectorable  thrusters  and  MSGCMG.  The  reference  path  for 
the  center  of  mass  of  the  simulator  consists  of  200  waypoints,  taken  at  angular  intervals 
of  1.8  deg  along  a  circle  of  diameter  lm  with  a  center  at  the  point  [2.0  m,  2.0  m]  in  the 
ICS,  which  can  be  assumed,  for  instance,  to  be  the  center  of  mass  of  the  target.  The 
reference  attitude  is  taken  to  be  zero  throughout  the  maneuver.  The  entire  maneuver  lasts 
147  s.  During  the  first  10  s,  the  simulator  is  maintained  fixed  in  order  to  allow  the 
attitude  Kalman  filter  time  to  converge  to  a  solution.  At  10  s  into  the  experiment,  the 
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solenoid  valve  regulating  the  air  flow  to  the  linear  air  bearings  is  opened  and  the 
simulator  begins  to  float  over  the  epoxy  floor.  At  this  point,  the  simulator  begins  to 
follow  the  closed  path  through  autonomous  control  of  the  two  thrusters  and  the 
MSGCMG. 

As  evidenced  in  Figure  7a  through  Figure  7d,  the  components  of  the  center  of 
mass  of  the  simulator  as  estimated  by  the  translation  linear  quadratic  estimator  are  kept 
close  to  the  reference  signals  by  the  action  of  the  vectorable  thrusters.  Specifically,  the 
mean  of  the  absolute  value  of  the  tracking  error  is  1 .3  cm  for  A rBI , ,  with  a  standard 

deviation  of  9.1  mm,  1.4  cm  mean  for  Arm  2  with  a  standard  deviation  of  8.6  mm,  2.4 
mm/s  mean  for  Av/;/ ,  with  a  standard  deviation  of  1.8  mm/s  and  3.0  mm/s  mean  for 
Avg 1 2  with  a  standard  deviation  of  2.7  mm/s.  Furthermore,  the  mean  of  the  absolute 
value  of  the  estimated  error  in  rBI  {  is  2  mm  with  a  standard  deviation  of  2  mm  and  4  mm 
in  rBI  2  with  a  standard  deviation  of  3  mm.  Likewise,  Figure  7e  and  Figure  7f  demonstrate 

the  accuracy  of  the  attitude  tracking  control  through  a  comparison  of  the  commanded  and 
actual  attitude  and  attitude  rate.  Specifically,  the  mean  of  the  absolute  value  of  tracking 
error  for  A6*3is  0.14  deg  with  a  standard  deviation  of  0.11  deg  and  0.14  deg/s  for 

A  ox  with  a  standard  deviation  of  0.15  deg/s.  These  control  accuracies  are  in  good 
agreement  with  the  set  parameters  of  the  Schmitt  triggers  and  the  LQR  design. 

Figure  8a  through  Figure  8d  report  the  command  signals  to  the  simulator’s 
thrusters  along  with  their  angular  positions.  The  commands  to  the  thrusters  demonstrate 
that  the  Schmitt  trigger  logic  successfully  avoids  chattering  behavior  and  the  feedback 
linearized  controller  is  able  to  determine  the  requisite  thruster  angles.  Figure  8e  and 
Figure  8f  show  the  gimbal  position  of  the  miniature  single-gimbaled  control  moment 
gyro  and  the  delivered  torque.  Of  note,  the  control  system  is  able  to  autonomously 
maneuver  the  simulator  without  saturating  the  MSGCMG. 


48 


Transversal  CoM  Position 


0.025 


Transversal  CoM  Velocity 


b) 


0.02 

0.015 

0.01 

0.005 

0 

-0.005 

-0.01 

-0.015 

-0.02 

-0.025 


0  20  40  60  80  100  120  140  160 

t  (sec) 


Longitudinal  CoM  Position 


Longitudinal  CoM  Velocity 


d) 


Z-axis  Attitude 


Z-axis  Attitude  Rate 


Figure  7.  Logged  data  versus  time  of  an  autonomous  proximity  maneuver  of  NPS 
SRL’s  3 -DoF  spacecraft  simulator  along  a  closed  path  using  vectorable  thrusters  and 
MSGCMG.  The  simulator  begins  floating  over  the  epoxy  floor  at  t  =  10  s. 
a)  Transversal  position  of  the  center  of  mass  of  the  simulator  in  I  b)  Transversal 
velocity  of  the  center  of  mass  of  the  simulator  in  ICS  c)  Longitudinal  position  of  the 
center  of  mass  of  the  simulator  d)  Longitudinal  velocity  of  the  center  of  mass  of  the 

simulator  e)  Attitude  f)  Attitude  rate 
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Figure  8.  Control  actuator  actions  during  autonomous  proximity  maneuver  of  NPS 
SRL’s  3-DoF  spacecraft  simulator  along  a  closed  path  using  vectorable  thrusters  and 
MSGCMG.  a)  Thruster  1  firing  profile  b)  Thruster  1  position  c)  Thruster  2  firing  profile 
d)  Thruster  2  position  f)  MSGCMG  torque  profile  e)  MSGCMG  gimbal  position 


Figure  9  depicts  a  bird’s-eye  view  of  the  spacecraft  simulator  motion.  Of 
particular  note,  the  good  control  accuracy  can  be  evaluated  by  the  closeness  of  the  actual 
ground-track  line  to  the  commanded  circular  trajectory  and  of  the  initial  configuration  of 
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the  simulator  to  the  final  one.  The  total  AV required  during  this  experimental  test  was 
0.294  m/s  which  correspond  to  a  total  impulse  of  7.65  Ns. 


1  1.2  1.4  1.6  1.8  2  2.2  2.4  2.6  2.8  3 

Xc(m) 

Figure  9.  Bird’s-eye  view  of  autonomous  proximity  maneuver  of  NPS  SRL’s  3-DoF 
spacecraft  simulator  along  a  closed  path  using  vectorable  thrusters  and  MSGCMG 

2.  Autonomous  Proximity  Maneuver  Using  Vectorable  Thrusters  Along 
a  Closed  Circular  Path 

Figure  10,  Figure  11,  and  Figure  12  report  the  results  of  maneuvering  the 
spacecraft  simulator  along  the  same  reference  maneuver  as  in  Section  6.1  but  by  using 
only  the  vectorable  thrusters.  This  maneuver  is  presented  to  demonstrate  the 
experimental  validation  of  the  STLC  analytical  results.  As  before,  during  the  first  10  s, 
the  simulator  is  not  floating  and  kept  stationary  while  the  attitude  Kalman  filter 
converges. 

The  tracking  and  estimation  errors  for  this  maneuver  are  as  follows  with  the 
logged  positions,  attitudes  and  velocities  shown  in  Figure  10.  The  mean  of  the  absolute 
value  of  the  tracking  error  is  1.4  cm  for  A rBI , ,  with  a  standard  deviation  of  8.5  mm,  1.4 
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cm  mean  for  ArBI  2  with  a  standard  deviation  of  8.6  mm,  2.5  mm/s  mean  for  AvBI ,  with  a 
standard  deviation  of  1.9  mm/s  and  3.1  mm/s  mean  for  A rBI2  with  a  standard  deviation  of 
2.8  mm/s.  The  mean  of  the  absolute  value  of  the  estimated  error  in  rBI ,  is  3  mm  with  a 
standard  deviation  of  3  mm  and  4  mm  in  rB1 2  with  a  standard  deviation  of  5  mm.  The 
mean  of  the  absolute  value  of  tracking  error  for  A 6i  is  0.52  deg  with  a  standard  deviation 
of  0.3 1  deg  and  0.24  deg/s  for  A with  a  standard  deviation  of  0.20  deg/s.  These  control 

accuracies  are  in  good  agreement  with  the  set  parameters  of  the  Schmitt  triggers  and  LQR 
design. 

Figure  11  reports  the  command  signals  to  the  simulator’s  thrusters  with  the 
commands  to  the  thrusters  again  demonstrating  that  the  feedback  linearized  controller  is 
able  to  detennine  the  requisite  thruster  angles  to  take  advantage  of  this  fully  minimized 
actuation  system.  Figure  12  depicts  a  bird’s-eye  view  of  the  motion  of  the  simulator 
during  this  maneuver.  The  total  AV  required  during  this  experimental  test  was  0.327  m/s 
which  correspond  to  a  total  impulse  of  8.55  Ns. 
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Figure  10.  Logged  data  versus  time  of  an  autonomous  proximity  maneuver  of  NPS 
SRL’s  3-DoF  spacecraft  simulator  along  a  closed  path  using  only  vectorable  thrusters. 
The  simulator  begins  floating  over  the  epoxy  floor  at  t  =  10  s.  a)  Transversal  position  of 
the  center  of  mass  of  the  simulator  in  the  I  b)  Transversal  velocity  of  the  center  of 
mass  of  the  simulator  in  ICS  c)  Longitudinal  position  of  the  center  of  mass  of  the 
simulator  d)  Longitudinal  velocity  of  the  center  of  mass  of  the  simulator  e)  Attitude 

f)  Attitude  rate 
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Thruster  1  Angle 


Figure  11.  Control  actuator  actions  during  autonomous  proximity  maneuver  of  NPS 
SRL’s  3-DoF  spacecraft  simulator  along  a  closed  path  using  only  vectorable  thrusters,  a) 
Thruster  1  firing  profile  b)  Thruster  1  position  c)  Thruster  2  firing  profile  d)  Thruster  2 

position 
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Figure  12.  Autonomous  proximity  maneuver  of  NPS  SRL’s  3-DoF  spacecraft 
simulator  along  a  closed  path  using  only  thrusters 
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III.  VECTORABLE  THRUSTER  CONFIGURATION  FOR  ROTO- 
TRANSLATIONAL  CONTROL  OF  SPACECRAFT 

A  novel  configuration  of  two  oppositely-mounted  vectorable  thrusters  is 
presented  to  provide  six  DoF  roto-translational  control  of  a  spacecraft.  Figure  13  depicts 
a  potential  vectorable  thruster  solution  based  on  the  research  presented  by  Canfield  and 
Reinholz  (1998).  By  allowing  each  thruster  to  be  vectored  within  a  hemispherical  space, 
small-time  local  controllability  of  underactuated  spacecraft  is  achieved.  This 
configuration  of  control  actuators  is  proposed  as  a  possible  alternative  to  conventional 
fixed  thruster-based  architectures  on  rapidly  deployable,  low-cost,  and  low-mass 
spacecraft  systems.  It  may  also  be  used  for  more  traditional  spacecraft  such  as  NASA’s 
Crew  Exploration  Vehicle  (CEV)  and  ESA’s  Automated  Transfer  Vehicle  (ATV).  The 
CEV’s  fixed  thruster-based  Reaction  Control  System  (RCS)  consists  of  a  configuration  of 
four  pods  of  six  thrusters  to  provide  six  DoF  roto-translation  control  (Jackson  and 
Gonzalez  2007).  In  comparison,  the  ATV  achieves  six  DoF  control  via  32  fixed  thrusters 
with  28  dedicated  to  the  Attitude  Control  System  (ACS)  and  4  dedicated  to  the  Orbital 
Control  System  (OCS)  (Cavrois,  Reynaud,  Personne,  Chavy  and  Strandmoe  2008).  The 
potential  advantages  of  reducing  the  number  of  thrusters  required  for  6-DoF  roto- 
translational  control  include  simplification  of  the  overall  spacecraft  design,  reduction  of 
the  required  fuel  for  a  desired  maneuver,  and  a  reduction  in  the  overall  size  of  the 
propulsion  system. 


57 


Figure  13.  Depiction  of  potential  hemispherical  vectorable  thruster  (Courtesy  Dr. 

Steven  Canfield,  Tennesse  Tech  University) 

This  chapter  presents  the  dynamics  model  of  this  novel  unconventional  spacecraft 
architecture  within  a  proximity  operations  environment  and,  following  the  mathematical 
methods  presented  in  the  previous  chapter,  the  small-time  local  controllability  of  this 
underactuated  system  is  demonstrated.  By  adding  two  additional  pairs  of  fixed  thrusters, 
the  system  becomes  fully  actuated  and  thus  input-output  linearizable.  Given  this  input- 
output  linearizability,  a  feedback  linearized  control  law  is  derived  to  that  can  control  the 
6-DoF  roto-translation  problem  via  traditional  linear  control  methods.  The  reference 
spacecraft  is  assumed  to  be  rigid  with  principal  moments  of  inertia  aligned  with  the  body- 
fixed  control  axes. 

A.  SPACECRAFT  RELATIVE  MOTION  DYNAMICS 

Consider  two  close-orbiting  spacecraft  moving  in  the  gravitational  field  of  Earth 
or  one  chaser  spacecraft  moving  with  respect  to  a  reference  point.  It  is  assumed  that  the 
chaser-target  pair  are  homogeneous  rigid  bodies  and  that  the  chaser-target  pair  are 
moving  in  a  nearly  circular  orbit  with  the  final  desired  state  of  the  control  problem  being 
a  soft  dock  of  the  two  vehicles. 
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Three  reference  frames  are  established:  an  inertial  frame  I  with  orthogonal  axes 
defined  by  the  set  of  unit  vectors  j  /' ,  i2 ,  /,  j ,  a  rotating  satellite  reference  frame  Jf  with 

orthogonal  axes  defined  by  the  set  of  unit  vectors  j  /?, ,  h2 ,  /?,  j  and  a  body-fixed 

framed?  with  orthogonal  axes  defined  by  the  set  of  unit  vectors  j bx , b2 ,  fx  }  .  The  rotating 

satellite  reference  frame  J-f  is  used  to  define  the  reference  coordinate  system  for  the 
chaser-target  rendezvous  problem  where  the  hx  is  collinear  with  the  position  vector  of  the 

target  spacecraft  with  respect  to  the  Earth’s  center,  h2  is  in  the  direction  of  the  velocity 
vector  of  the  target  spacecraft  aligned  with  the  local  horizontal  and  h}  is  normal  to  the 
orbit  plane.  In  this  manner,  the  angular  velocity  of  J-f  with  respect  to  I  is 

(a  hi  —  ~cofi2  (92) 

where  co0  is  the  target  spacecraft’s  mean  motion.  The  angular  velocity  of  f>  with  respect 
to  I  is  then  given  by 

(0BI  —  Q>bh  +  (oH1  —  —  co()h2  (93) 


1.  Translational  Motion  of  the  Chaser-Target  Pair 

The  3 -DoF  translational  motion  of  the  rigid  chaser  spacecraft  with  respect  to  3~f 
are  given  by  Vallado  (2001): 


rBH  ~'VBH 

V  BH  =  f  (  VBH  ’  V  BH  )  +  ^ 


(94) 


where  v(/)  e  Vv  =  ['-’«//, i,  y«//,2,  y/jw ,3]  is  the  translational  velocity  vector  of  the 
chaser’s  body-fixed  frame  J7  with  respect  to  the  moving  satellite  reference  frame  Jf 
expressed  in  ,  r  (7)  e  M3 ,  Vr  =  \jBH , ,  rBH  2 ,  rBH  3  ]  is  the  position  vector  of  the  chaser’s 
body-fixed  frame  2?  with  respect  to  the  moving  satellite  reference  frame  2>  expressed 
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in  22 ,  ffF(/)eR3,VFF  =  ^^,  V2,V3j  is  the  control  vector  acting  on  2f ,  m  is  the 
mass  of  the  chaser  spacecraft  and  f  (rBH ,  \  BH  )el3  represents  the  drift  vector 


f  {^BH’  VS//) 


^0J()rBH,\  ^■G)0VBH,2 

—2co0vbhi 

-co2r 

^0  r  BH, 3 


(95) 


2.  Rotational  Motion  of  the  Chaser  Spacecraft  with  Respect  to  the 
Inertial  Frame 

The  rotational  motion  of  the  chaser  spacecraft  about  body-fixed  axes  with  their 
origin  located  at  the  body’s  center  of  mass  can  be  described  using  Euler’s  equations.  For 
the  general  case,  in  which  the  control  axes  do  not  coincide  with  the  principal  axes  of 
inertia  and  torquing  devices  such  as  independent  gas  thrusters  are  used  to  impart  torques 
about  these  control  axes,  these  become  from  Wie  (1998,  p.  341): 

<aBI  =  -J-'ofBIJo)BI  +  J-1  T  (96) 

where  <oBI  (/  )  e  R3 ,  Vas/  =  \joBIX,coBI2,  coBI  3  ]  is  the  angular  velocity  vector  of  the  body- 
fixed  frame  22  with  respect  to  an  inertial  frame  I  expressed  in  22 , 
T (t)  e  R3,  VT  =  [  bT{,  bT2,  is  the  control  torque  vector  acting  on  22 .  The  notation 
<x>xBI  e  R3'3  denotes  the  skew-symmetric  matrix 
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is  the  inertia  matrix  with  respect  to  the  J2  given  by 
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3.  Quaternion  Kinematic  Representation  of  the  Chaser  Spacecraft’s 
Orientation 

The  orientation  of  a  rigid  spacecraft  can  be  defined  through  various 
parameterizations  of  the  special  orthogonal  group  50(3) .  In  this  work,  the  quaternion  is 

considered  to  parameterize  this  orientation  due  to  their  suitability  for  onboard  real-time 
computation  and  subsequent  common  use  for  spacecraft  attitude  determination  (Wie  et  al. 
1989).  The  unit  quaternion,  first  envisioned  by  Hamilton  in  1843,  can  be  separated  into 
two  parts,  the  first  being  a  three-dimensional  vector  indicating  the  direction  of  the  axis  of 
rotation  and  the  second  being  a  scalar  quantity  which  relates  the  angle  of  rotation  about 
this  axis  of  rotation  such  that 


q  =  q  +  q4=q1i  +  q2j  +  q3k  +  q4 
where  q(t)  =  |q(t),^4(t)}elR3vlR  and 


q{  =  bt  sin(J/2),  i  =  1,2,3 
qA  =cos(J/2) 


(99) 


(100) 


where  5  represents  the  magnitude  of  the  rotation  around  the  chosen  axis  and 
( bx ,  b2 ,  b3 )  arc  the  direction  cosines  locating  the  axis  of  rotation  to  the  inertial  frame  by 
Euler’s  Theorem  (Wie  and  Barba  1985;  Wie  et  al.  1989;  Cristi  and  Burl  1993).  The 
properties  of  the  unit  vectors  (i,  j,  k)with  respect  to  the  chosen  reference  frame  when 
multiplied  are  such  that  0 

ij  =  -ji  =  k,  jk  =  -kj  =  i,  ki  =  -ik  =  j,ii  =  jj  =  kk  =  -1  (101) 

By  considering  the  direction  cosines  expressed  in  (100)  to  be  equivalent  to  the 
unit  vectors  (i,  j,  k) ,  the  quaternion  can  be  used  to  define  a  sequence  of  rotations  about  a 
series  of  Euler  angles  defined  by  (j),6 ,  and  y/  .  By  considering  a  3-2-1  body-fixed 
rotation  sequence  from  2?  to  I ,  as  is  typical  for  visualizing  spacecraft  orientations,  a  set 
of  elementary  quaternions  can  be  defined  from  the  three  Euler  angles  similar  to  the  three 
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sets  of  elementary  rotation  matrices  (Wie  1998).  In  this  manner,  8X  =  <f>,S2  =  0,S3  =i//  and 
from  (100)  and  (101),  the  quaternion  associated  with  this  sequence  of  rotations  becomes 


9b/, i 

sxc2c3  -C,s2s3 

9b/, 2 

cxs2c3  +  s{c2s3 

9b/,  3 

Cl C2s3  -svs2c3 

<a> 

V 

1 _ 

_clc2c3+sls2s3_ 

(102) 


where  st  =  sin (cf/2), c(.  =  cos (<5)/2),i  =  1,2,3 .  Note  that  the  order  of  elementary 

quaternion  multiplications  is  reversed  from  the  order  of  elementary  rotation  matrices  due 
to  the  properties  of  quaternion  multiplication  expressed  in  (101). 

The  kinematic  differential  equation  of  the  quaternion  representing  the  orientation 
of  the  chaser  spacecraft’s  body-fixed  frame  %  to  7  is  given  as 


9.B/  2  (  9 B/, 4^3x3  9  b/  )  ®*B/ 

'  =-I  r 

9b/,4  2  BI 


(103) 


where  qxBI  e  R 3,3  represents  the  skew-symmetric  matrix 


9b/ 


0  “9b/,3  QbI,! 

^BIA  0  ~ 9b/,1 

_ 9b/,  2  9b/,1  0 


(104) 


and  f3x3  represents  a  3x3  identity  matrix.  One  can  quickly  show  using  (100)  that  the 
quaternion  is  subject  to  the  constraint 

9b/, 4  =  1 —  Qb/Qb/  =  1 —  |*Ib/ 1  (105) 


62 


4.  Angular  Motion  of  the  Chaser  Spacecraft  with  Respect  to  the  Satellite 
Reference  Frame 

Considering  the  rotating  satellite  reference  frame  3~f ,  the  quaternion  describing 
the  orientation  of  £  with  respect  to  Jf  is  defined  as  above  so  that 
qBH  (0  =  (7),  <7 bh  4  (0}  e  •  From  (103),  the  kinematic  differential  equation  of 

q bh  is  §iven  by 


q  BH  n.  (  QbH,  4^3x3  QbH  )(>iBH 


tfBHA  2  q  BH®*BH 


(106) 


Likewise,  the  constraint  on  the  components  of  qBH  hold  so  that  (105)  becomes 


y  bh  a  i  q^ 


The  3x3  rotation  matrix  that  brings  J~f  onto  2>  is  given  by 


R h  [qBH a  qJI  +  2qs//qB//  2qBHAqBH 


so  that 


1 

_ 1 

1 

_ 1 

- 

k 

=  hRb 

4 

= 

i 

i _ 

i _ 

- 

^11  ^21  ^31 

Rn  R22  Rn 

R\i  ^23  ^33 


(107) 


(108) 


(109) 


From  (93)  and  (108),  the  angular  velocity  of  the  body-fixed  frame  2?  with 
respect  to  J~f  expressed  in  2>  is  therefore 


^ HI  ®8/  [0,0, 

where  (t) 6  R  ,\/o)BH  —  \j^,bha,^>bh,2,^>bha\ 


(110) 


Taking  the  time  derivative  of  (1 10)  yields 
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(Ill) 


which,  given 


where  ofBH  e  M3*3 


®  BH  ^  BI  [Oj  0,  0\  ] 


si?„[0,0,n>0f  =  ofBHRH  [0,0,  <y„]7 


represents  the  skew-symmetric  matrix 


®BH 


0  ~^BH,  3  (°BH,1 

®BH,3  0  ~(°BH,  1 

~G)BH,  2  ^RH  ,\  0 


yields 


to 


—  6)  m  +  G)bhRh 


[0,0,  (tf0]7 


(112) 


(113) 


(114) 


Therefore,  by  substituting  (96)  into  (114),  the  angular  motion  equations  of  the 
chaser  spacecraft’s  body-fixed  frame  2> with  respect  to  the  rotating  satellite  reference 
frame  3~f  expressed  in  2>  become 


*\bH  2  \^BH  ,4^3x3  +  ^BHj^BH 

1  T 

^BHA  =  _  (115) 

6)bh——J  (-'J/j/  /a  ~  (’hfihiH  [  3  ’  ^2i  ’  ^33  ]  +  J  T 


5.  Chaser  Spacecraft  Relative  Motion  Dynamics 

From  (94)  and  (115),  the  full  set  of  equations  describing  the  relative  motion  of  the 
chaser  spacecraft  with  respect  to  the  rotating  satellite  reference  frame  Jf  become 
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rBH  V  BH 


^  BH  f  (^BH  ’  V«W  )  +  fYl  Rb F 
O-RH  =  ^{‘IbHA^xI  +  ^BH^^BH 


(116) 


Qbha  2  ^BH^BH 

®  BH  —  —J  ^Bl^^BI  _  ^O^BH  [^13 ’^23 ’^33]  ^  ^ 

with  F  (7)  e  M3 ,  VF  =  ^ 11  Ft ,  BF2 ,  BF3  J  representing  the  control  thrust  vector  acting  on  the 
body-fixed  frame  2>  and  " RB  =  BRTH  representing  the  3x3  rotation  matrix  that  brings 
2?  onto  J-f . 

Following  Wie  (1998,  p.  407),  the  system  of  concern  can  be  reduced  to  the  twelve 
dimensional  space  through  substitution  of  the  constraint  on  the  components  of  the 
quaternion  given  by  (105)  into  (108)  and  (1 16)  to  yield 


rBH  V BH 


V BH  ~  f  (  rBH  ’  V  RU  )  +  m  RB^ 


*\bH  2  Ifis/zl  1 1x1  +  QbH  j  ^BH 

6)bh——J  cogiJ(oBI  —  co0coBH  [f?13,f?23,f?33]  +  J  T 


(117) 


where 


Brh  =  (i  “  ■ 2  ||q«r  If )  hxi  +  2q  bh^bh  +  2^~\K 


Qbh 


(118) 


By  considering  a  diagonal  zed  inertia  matrix  so  that  (98)  simplifies  to 


J  = 


Ju  0  0 

0  J22  0 

0  0  J33 


(119) 


the  relative  motion  equations  described  by  (117)  can  be  compactly  expressed  in  proper 
control-affine  form  as 
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x  =  f  (x)  +  G(x)u 

where  x(t)  e  R12  is  the  12-dimensional  state  vector  given  by 


(120) 


X  [x,,X2,X3,...,X12]  \rBH,l’rBH,2’rBH,3  ’  QbH,  1  2  ’  ^BH, 3  ’  VBH,  1  ’  ’  VBH ,3  ’  ^BH,2  ’  ®S/f , 3  ]  (121) 

u(0  e  U6  is  the  6-dimensional  control  vector  given  by 


[n1,n2,...,n6]r  =  BFl,BF2,BFi,  %,BT2,% 


(122) 


the  drift  vector  f  (x)  e  R12  given  by 


f(x)  = 


-  —  x  [x  +  co  (2 rjx  +  2x  x  )]  +  —  x  [x  +  co  (l  -  2x '  -  2x"  )1  +  —  q  [x  +®  (2xx  -  2 rjx  )] 
2  °  4  2  4  2 

~X  [-T,0  +  (2X4X,  “  2/7X5  )]  +  _X4  [XP  +  CO0  (i  -  2X4  “  2xi  )]  [X„  +  (2?7X4  +  2X5X6  )] 

2  2  2 

-  —  x5  [x„  +  a>Q  (2x4x6  -  2 rjxs )]  +  —  x4  [xn  +  coa  (2 r/xt  +  2xsx6 )]  +  —  rj  [xp  +  co..  ( 1  -  2x~  -  2x4  )] 


3  co  x  +2  co  x 

0  1  0  8 


O',,"',,) 


O',,-',,) 


(j  -J  ) 

\  B,  1  B ,2  / 


[xp  +  ®o  (l  -  2xJ  -  2x(  [x,  +  <»  (2t;x4  +  2xsx6  )]  +  <»  [xn  (l  -  2x"  -  2x_  )  -  xp  (2  r/xt  +  2xsxs )] 

[x,2  +  £9  (I  -  2x)  -  2x‘ )]  [x„  +  £9  (2x4X6  -  2t7X5  )]  +  «  [xp  ( 2x4X6  -  2 jjxs  )  -  xio  (l  -  2x~  -  2x‘ )] 
[x,  +  £9o  (2  r/xt  +  2xxs  )][x„  +  a>o  (2x4X6  -  2r/x j )]  +  a>o  [x  (2t;x4  +  2xixs )  -xn  (2x4X6  -  2  r/xs )] 

(123) 


where  77  =  ±Jl-  qBH  and  the  control  matrix  G(x)  e  R12*6  given  by 
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(124) 


G(x) 


where  Gx  (x)  e  R6'6  represents  the  control  distribution  matrix  given  the  placement  of  the 

control  inputs  with  respect  to  the  spacecraft’s  mass  center.  For  the  ideal  case,  where  each 
thrust  vector  is  parallel  to  its  corresponding  body-fixed  axis  acting  through  the 
spacecraft’s  mass  center  and  three  independent  paired  thrusters  provide  torque  about  their 
corresponding  body-fixed  axes,  Gj(x)  becomes 


“  RB  °3.v3 

03,3  J  X  _ 


(125) 


B.  SMALL-TIME  LOCAL  CONTROLLABILITY  CONSIDERATIONS  FOR 
THE  PROPOSED  6-DOF  SPACECRAFT  DESIGN 

Small-time  local  controllability  for  the  6-DoF  nonlinear  equations  representing 
the  relative  motion  of  a  chaser  spacecraft  given  by  (120)  can  be  determined  following  the 
same  methods  as  presented  in  Chapter  II  for  the  3 -DoF  spacecraft  simulator.  In  order  to 
properly  deal  with  the  momentum  of  the  drift  vector  term  f(x),  the  analysis  of  small¬ 
time  local  controllability  is  performed  at  the  equilibrium  point  corresponding  to  f  (p)  =  0 
or 


P 


[0,  x2 , 0, 0, 0, 0, 0, 0, 0,  o,  o,  o]7 


(126) 


By  varying  the  control  vector  u ,  the  minimum  number  of  control  actuators  for  the 
system  can  be  gained.  In  order  to  consider  a  general  case,  the  thrust  vector  in  the  b] 
direction  is  defined  as  BFX  and  is  taken  to  be  directed  through  the  spacecraft’s  mass 
center  so  that  no  torque  is  imparted  on  the  body;  the  thrust  vector  parallel  to  the  b2  axis  is 
defined  as  BF \  rot3  and  is  taken  to  be  a  distance  L  from  the  spacecraft’s  mass  center  along 
the  bx  axis  so  that  both  a  force  is  imparted  in  the  b2  direction  as  well  as  a  torque  about 
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the  b3  axis;  and  the  thrust  vector  parallel  to  the  b2  axis  is  defined  as  BF;  rot2  and  is  also 
taken  to  be  a  distance  L  from  the  spacecraft’s  mass  center  along  the  bx  axis  so  that  both 
a  force  is  imparted  in  the  b3  direction  as  well  as  a  torque  about  the  b2  axis.  Additionally, 
the  effects  of  pure  control  torques  are  considered.  Each  torque  is  assumed  to  provide  a 
pure  torque  about  its  corresponding  body-fixed  axis  so  that  HT]  represents  a  torque 

about  hj,  bT2  represents  a  torque  about  b2,  and  BT3  represents  a  torque  about  b2 .  In  this 
manner,  a  small-time  local  controllability  analysis  is  conducted  on  variations  of  the 
control  vector  u(t)  e  U6,  Vu  =  \_BFl,BF2/rot3,  BF3/rot2,  BTl,  BT2,  5T3]  such  that  the 

corresponding  vectors  of  the  associated  control  distribution  matrix  G1  (x)  are: 


bF]  ->  g (x)  =  [ m1Rn , m  lRn, m  lRn , 0, 0, 0]r 

(127) 

BF2,ro,3  ->  g  (X)  =  [m~lR21  > 111  'R22 ,  m  '  R23  >  0- 0,  Ut  ]F 

(128) 

BFyrot2  ->  g (x)  =  \_mxRn, m  'R,2 , m  ' R,, , 0, LJ2  , oj 

(129) 

^^g(x)  =  [o,o,o,zjr1,o,o]r 

(130) 

i,r2->g(x)  =  [o,o,o,o,zj2-I,o]r 

(131) 

%  ->  g  (x)  =  [0,0,0, 0,0,LJ?J 

(132) 

The  results  of  this  investigation,  which  are  summarized  in  Table  5,  demonstrate 
that  the  minimum  number  of  actuators  for  small-time  local  controllability  of  the  relative 
motion  spacecraft  dynamics  defined  by  (1 17)  include  one  thrust  vector  acting  through  the 
spacecraft’s  center  of  mass  and  two  body-fixed  control  torques.  This  combination  can  be 
achieved  by  considering  a  total  of  six  control  actuators  with  two  fixed  thrusters  mounted 
on  opposite  faces  of  the  spacecraft  with  their  thrust  vectors  acting  in  an  opposing  manner 

along  the  b]  axis  through  the  spacecraft’s  center  of  mass  and  two  sets  of  fixed  paired 
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thrusters  acting  together  to  produce  pure  torque  about  two  of  the  body-fixed  axes.  For 
instance,  by  considering  the  control  vector  u(/)eU  3,Vu  =  [BFl,BT2,BT3~\  and  the 

corresponding  control  distribution  matrix  Gx  (x)  e  M6a3 


G1(x)  =  [g1,g2,g3] 


mlRn  0  0 

m'Rn  0  0 

mlRu  0  0 

0  0  0 

0  LJ~l  0 

0  0 


(133) 


the  Lie  Algebra  evaluates  to: 


./’(A)  =  span  < 


g,,g2,g3,[f,g1],[f,g2],[f,g3],[f,[f,g1]],[f,[f,g2]], 

[go[f5g2]]>[go[f5g3]]{g25[f>g3]][f{t[f,gi]]] 


(134) 


However,  by  considering  the  vectorable  thruster  depicted  in  Figure  13,  only  two 
control  actuators  are  required  to  achieve  small-time  local  controllability  of  the  relative 
motion  spacecraft  dynamics  given  by  (117).  This  can  be  achieved  by  placing  each 
hemispherically  vectorable  thruster  a  distance  L  from  the  spacecraft’s  center  of  mass 

along  the  bx  axis.  For  this  configuration,  the  control  vector 

u(t)  e  U3,  Vu  =  [sFi,  BF2/rot3,  BFyrotl J  is  considered  with  the  corresponding  control 
distribution  matrix  G,  (x)  e  IR6*3 


GI(x)  =  [g1,g2,g3]  = 


m 

'Rn 

m  R2[ 

m  «3I 

III 

%2 

mlR22 

m  'Rn 

III 

X 

m  lR2i 

m  lR33 

0 

0 

0 

0 

0 

LJ-] 

0 

LJ~l 

^33 

0 

(135) 


yielding  a  Lie  Algebra  for  this  minimally  configured  system  of 


69 


£  (A)  =  span 


gl>g2?g3>[f5g1]5[f5g2]5[f>g3]>[t[f,g1]],[f,[f,g2]], 

[  g  1  ^  [  f  ’  g  3  ]  ]  ^  [  g  2  ^  [  f  ’  g  3]],[f,[  g  1  ^  [  f  ’  g  3  ]]]  9  [f’  [  g  2  9  [  ’  g  3  ]  ]  ^ 


(136) 


Control  Inputs 

dim(£(A)) 

Controllability 

bF  or  bF  or  BF 

1  1  ur  1  2/ rot 3  Ur  1  3/rot2 

3 

Inaccessible 

(AAm.jMAA.M 

6 

Inaccessible 

(AC 

8 

Inaccessible 

( “F, ,  T2  ) ,  ( ‘F, ,  T,  ) ,  ( A* .  T,  ) ,  ( A™, .  T,  ) , 
('F>M2.Tl  )"'(  A,,,...7"'  ) 

9 

Inaccessible 

( BF2,ro,3 ,  BF3lrot2  ) ,  ( BF2/rot 3 ,  Tx  )  or  ( BF3/rot2 ,  Tx  ) 

10 

Inaccessible 

(BFx,BTx,BT2),(BFx,BTx,BT3)or(BFx,BT2,BT3) 

12 

STLC 

B  Z7  B  tj  B  tj 

r\  9  r 2! rot2  9  r2!rot2 

12 

STLC 

Table  5.  STLC  Analysis  for  the  6-DoF  Relative  Motion  Spacecraft  Dynamics 


C.  PROPOSED  FEEDBACK  LINEARIZABLE  SPACECRAFT  SYSTEM 

DESIGN  WITH  MINIMUM  NUMBER  OF  CONTROL  ACTUATORS 

In  the  previous  section,  a  full  analysis  of  the  small-time  controllability  for  a 
spacecraft  with  variations  on  the  control  inputs  was  conducted.  This  analysis  yielded  two 
satisfactory  variations  in  control  inputs  for  small-time  local  controllability:  a  set  of  six 
control  actuators  consisting  of  two  opposing  fixed  thrusters  acting  through  the 
spacecraft’s  center  aligned  with  one  of  the  three  body-fixed  axes  and  two  sets  of  paired 
thrusters  acting  about  two  of  the  three  body-fixed  axes  or  a  set  of  two  control  actuators 
consisting  of  oppositely  mounted  hemispherically  vectorable  thrusters.  From  this 
analysis,  a  unique  spacecraft  control  actuator  configuration  is  proposed  that  leverages  this 
study  to  provide  a  feedback  linearizable  system  as  might  be  necessary  in  a  path  or  sensor 
constrained  environment.  As  stated  in  Section  II,  in  order  for  a  system  to  be  state- 
feedback  linearizable,  it  must  be  square  so  that  the  number  of  actuators  is  equal  to  the 
number  of  observed  states.  From  this,  the  number  of  observed  states  and  number  of 
control  inputs  are  assumed  to  be  six  and  a  novel  minimally  actuated  system  is  proposed 
consisting  of  a  total  of  four  thrusters.  Figure  14  graphically  depicts  the  forces  acting 
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upon  the  body-fixed  frame  2?  through  two  hemispherically  vectorable  oppositely 
mounted  thrust  vectors  Fa  e  M3  and  F,,  e  E'  and  a  set  of  paired  thrusters  acting  about  the 

bx  axis  to  impart  the  torque  fe  1 . 


k 


Figure  14.  6-DoF  Minimally  Actuated  Feedback  Linearizable  Spacecraft 

Configuration 

In  order  to  properly  resolve  Fa  and  FA  into  2? ,  the  vectorable  thruster  frame 
Tfis  established  with  orthogonal  axes  defined  by  the  set  of  unit  vectors  { vn ,  vi2 ,  vi3 } 
where  i  =  a,b  corresponding  to  either  Fa  or  F^ .  Both  axes  va2  and  va2  are  taken  to  be 
collinear  with  the  b1  axis.  Additionally,  the  va3  axis  is  taken  to  be  aligned  with  the  bx 

axis  with  va3  =  bt  while  the  vM  axis  is  taken  to  be  aligned  with  bt  with  v/)3  =  -b{ .  The 
coordinates  of  the  imparted  force  vector,  as  depicted  in  Figure  14,  is  defined  in  the 
thruster  frame  through  two  angles  «.  and  /? ,  referred  to  as  the  polar  and  azimuthal 

angles,  respectively.  The  azimuthal  angle  ai  is  the  angle  between  the  axis  of  the 
imparted  thrust  vector  and  the  v.3  unit  vector  and  is  limited  to  the  range  [0,^/2] .  The 
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polar  angle  is  defined  as  the  angle  between  the  vn  and  the  projection  of  the  thrust  vector 
F(  onto  the  {vn,vi2}  plane.  In  this  manner, 


and 


'  Fa  =  ~Fa  [sin  aa  cos  fia ,  sin  aa  sin  /3a ,  cos  a„  f 


'  Fh  =  -Fb  [sin  ah  cos  J3h ,  sin  ah  sin  (5h ,  cos  ah  f 


The  3x3  rotation  matrices  that  bring  IF  onto  2?  are  then 


(137) 


(138) 


0  0  1 

0  1  0 

-10  0 


(139) 


and 


0  0 
0  1 
1  0 


-1 

0 

0 


(140) 


In  this  manner,  the  hemispherically  vectorable  thrust  vectors  can  be  expressed  in 
the  body-fixed  frame  2>  as 

B^a  -  ~Fa  [cosaa,sinaa  sin/?a,-sinaa  cos/?a]r  (141) 

and 

=  -Fb  [-cosaA,sinar6  sin  /^^in  ah  cos  J3hf  (142) 

yielding  the  thrust  vector  given  in  (1 17)  by  F (7)  e  R3,  VF  =  BFa  +  sFfc . 

The  torque  vector  in  (1 17)  given  by  Tel3  can  be  resolved  as 
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(143) 


-Fa  sin  au  cos  pa  -  Fh  sin  ah  cos  pb 
-Fa  sin  cta  sin  Pa  +  Fh  sin  ab  sin  Ph 


Defining  the 

control 

inputs 

as 

such 

yields 

the 

control 

vector 

u(/)  e  U‘,  Vu  =  [  ‘FM2,  X  • 

T2 ,  bF  ]  and  corresponding  control  distribution 

matrix  Gt  (x)  e  Mfiv6 

m  lRn 

zzz  R2l 

m  R, , 

0 

m  1 R, , 

zzz  xR2X 

m'Rn 

mlR22 

m%2 

0 

ZZZ  'i?32 

m  lR22 

G,M  = 

m  '  R]3 

0 

m1R23 

0 

mlR 13 

0 

0 

dJ\\ 

m'  R33 

0 

m  R23 

0 

(144) 

0 

0 

dd 22 

0 

—dJ2  2 

0 

0 

-dJ~l 

0 

0 

0 

dJ33 

The  components  ^  represent  the  entries  in  the  corresponding  z'-th  row  and  j- th 
column  of  11 RB  as  defined  by  (109)  and  the  variable  d  is  defined  similar  to  what  is 
proposed  for  the  3-DoF  spacecraft  simulator  presented  in  Chapter  II  whereby  the  control 
variables  are  defined  given  a  desired  control  input  to  the  system  with  respect  to  p  .  Thus 
if 


U  <  o  ( 


tiF, 


BU,  >  0 


HF 


d  =  -L 


~Fa  cos  ,  BF2lrot 3  =  -Fa  sin  sin  pa ,  BFVrot2  =  Fa  sin  cos  pa 
X  =  Tc ,  bT2  =  -Fb  sin  ab  cos  ph ,  BT2  =  -Fb  sin  ab  sin  ph 

d  =  L 


(145) 


Fh  cos  ab ,  BF2/rol 3  =  -Fb  sin  ab  sin  ph ,  BFVrot2  =  -Fb  sin  ah  cos  pb 
X  =  ~Tc  »  BTl  =  Fa  Sin  a  a  C0S  Pa  >  X  =  ~F a  Sin  Sin  Pa 


1.  Smooth  Feedback  Control  via  State  Feedback  Linearization 

The  methodology  presented  in  (62)  through  (65)  can  be  used  to  demonstrate  that 
the  nonlinear  relative  motion  dynamics  given  by  (117)  are  state  feedback  linearizable. 
With  the  nonlinear  relative  motion  dynamics  expressed  in  control-affine  form  given  by 
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(120)  with  Gj  (x)  defined  by  (144),  the  matrix  E (x)  as  defined  by  (63)  is  invertible  and 
the  relative  degree  of  the  system  at  the  equilibrium  point  p  is 
(rj,r2,r3,r4,r5,r6)  =  (2, 2, 2, 2, 2, 2) .  Therefore,  the  total  relative  degree  of  the  system  is 

twelve.  Given  that  the  total  relative  degree  of  the  system  is  equal  to  the  number  of  states, 
the  nonlinear  system  can  be  fully  linearized  by  state  feedback  to  yield  smooth  feedback 
control.  If  we  choose 

h(x)  =  [h!  (x)’h2  (x)]r  =[[xi,X2,x3],[x4,X5,x6]]r  =[rBH,qmf  (146) 
the  state  transformation  can  be  chosen  as 


That  is 


z  =  [hj  (x),h2  (x),Zfhj  (x),Zfh2  (x)J 


1, 

/ 

/  ||  m2  _  ' 

^BH  ’  QbH  ’  V BH  ’  ^  I 

\ 

/I  |  QjS//  ||  ^3x3+(1bH 

\®*BH 

(147) 


rBH  ~  Z1 
*\.BH  =  Z2 

vbh=  z3  (148) 

01 BH  —Q  Z4 

©a,  =  Q~\  +  BRH[0,0,a>0f 

where  z,  e  eR\z3eM3,z4eR3  are  new  state  variables  and  Q  e  M3*3  is  given  by 


0 


2 


(149) 


and  the  system  is  transformed  into 
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(150) 


Z1  =  Z3 


z3  =f(zi,z3)  +  //i?i;F 

\j-ls(Q-\  +  bRh  [0,0,fflbf )  j(q\  +  bRh  [0,0,®/) 
[5(e"1z4)%[0,0,fflb]r+J-1T 


i4=QQ'z4  +  Q 


+ 


where  5’  ( • )  represents  the  skew-symmetric  matrix  of  •  . 

The  relative  motion  dynamics  given  by  (117)  can  now  be  transformed  via  (150) 
and  the  state  feedback  control  law 

[F,T]r  =£'(x)~1(v-b)  (151) 


into  the  linear  system 


^6*6 

^ 6x6 

z  + 

^6x6 

_06.V6 

o6x6: 

/ 

L  6x6  J 

(152) 


where  fi  (x)  e  R6*6  =  ZGLfh(x),b  e  R6  =  Zfh(x)  and  v(t)  e  R6,  Vv  =  [v15  v2]r. 

With  the  system  now  in  linear  form,  traditional  linear  control  techniques  such  as 
the  linear  quadratic  regulator  as  presented  in  Chapter  II  can  be  used  to  stabilize  it  to  a 
desired  state. 


2.  Determination  of  the  Commanded  Azimuthal  and  Polar  Thruster 
Angles 

Using  the  decoupling  law  as  expressed  in  (151),  the  commanded  azimuthal  and 
polar  angles  for  the  hemispherically  vectorable  thrusters  and  the  required  torque  for  the 
paired  thrusters  can  be  determined.  For  example,  if  BUX  >  0  with  the  azimuthal  angle 
0  <  ab  <  nj 2  and  the  polar  angle  0  <  <  2n ,  the  angles  to  slew  thruster  b  can  be 

resolved  as 
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(153) 

(154) 

(155) 

(156) 

(157) 

(158) 


Likewise,  for  the  case  where  BU \  <  0  with  the  azimuthal  angle  0  <  ah  <  ir/l  and  the  polar 
angle  0  <  f5h  <  In ,  the  angles  to  slew  thruster  a  can  be  resolved  as 


Pb=atm2(BF2/rot3,-BF3lrot2) 


and  the  required  thrust  as 


Bp2  +  bf2  _i_  Bp?2 

1  1  ^  1  2/rot3  ^  1  Vrotl 


The  angles  to  slew  thruster  b  can  be  resolved  as 


(159) 

(160) 

(161) 
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(162) 


n 


Pa  =  atan2  (  BT3 ,  BT2 ) 


(163) 


and  the  required  thrust  as 


B  77  /  Brji 2  I  Brpl 


(164) 
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IV.  QUATERNION  FEEDBACK  ATTITUDE  REGULATION  OF 
UNDERACTUATED  RIGID  SPACECRAFT 

Methods  of  providing  three-axis  stabilization  via  two  body-fixed  control  torques 
are  introduced  in  this  chapter.  The  current  state-of-the-art  in  three-axis  stabilization  for  a 
rigid  spacecraft  with  arbitrary  inertia  via  two  body-fixed  control  torques  is  limited  to  rate 
stabilization  followed  by  three  sequential  Euler  angle  rotations  or  time -varying  attitude 
control  based  on  a  variety  of  switching  logics.  Byrnes  and  Isidori  (1991)  demonstrated 
that  for  the  full  underactuated  rigid  spacecraft  system  of  angular  motion  equations,  a 
smooth,  time-invariant  state  feedback  control  can  be  designed  to  locally  stabilize  the 
spacecraft  with  a  revolute  motion  about  the  unactuated  axis.  Without  violating  their 
assertion,  this  chapter  introduces  a  novel,  time-invariant  smooth  state  feedback  controller 
based  on  generalized  inverse  methods  and  quaternion  kinematics  to  obtain  simultaneous 
attitude  stabilization  to  an  arbitrarily  small  region  of  the  zero  error  state  with  two 
bounded  body-fixed  control  torques  for  an  asymmetric  spacecraft.  This  smooth  state 
feedback  controller  will  be  shown  to  be  most  naturally  suited  to  either  rest  to  rest 
reorientation  maneuvers  or  detumbling  and  reorientation  maneuvers.  It  is  not  applicable 
for  attitude  mainentenance  for  attitude  tracking  or  in  the  presence  of  large  disturbance 
torques  as  a  perturbation  is  imparted  about  the  unactuated  axis  through  the  interactions 
between  the  two  actuated  axes  to  stabilize  the  attitude  to  a  desired  orientation.  This 
chapter  provides  two  key  contributions:  First,  the  method  of  generalized  inverse  control 
is  extended  from  a  study  of  only  the  dynamics  as  considered  by  Bajodah  (2009b)  to  the 
full  system  of  dynamic  and  kinematic  equations.  Second,  two  separate  null-control 
vector  designs  are  derived  for  the  proposed  quaternion  feedback  regulator  which 
capitalize  on  the  generalized  inverse  control  methodology.  The  first  design  is  Lyapunov 
function  based  and  is  conjectured  to  yield  global  stability  but  results  in  significant  control 
effort  noise  near  an  inherent  singularity  in  the  null-control  vector.  The  second  design  is 
based  on  perturbed  feedback  linearization  and  provides  local  stability  while  keeping  the 
control  effort  to  low  levels.  In  addition  to  formal  proofs  for  each  control  design, 
simulation  results  are  provided  to  demonstrate  the  capabilities  of  these  novel  controllers. 
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A.  UNDERACTUATED  RIGID  SPACECRAFT  REORIENTATION 
THROUGH  QUATERNION  FEEDBACK 

In  this  section,  the  general  case  of  an  underactuated  rigid  spacecraft  being  rotated 
by  two  body-fixed  control  torques  is  considered.  In  order  to  simplify  the  problem,  an 
ideal  control  torque  is  assumed.  The  rotational  motion  and  quaternion  kinematic 
representation  of  the  orientation  of  the  rigid  body  remain  as  developed  in  Section  III  for 
the  general  case  of  the  body-fixed  frame  %  with  respect  to  an  inertial  frame  I .  In 
summary,  these  are: 


d)B1  =  -J  1(D*1J(dm  +  J  'T 

Qm  =  QbI  (165) 

•  __i  T 

VbIA  2 

1.  Attitude  Error  Parameterization  via  Quaternions 

The  3x3  rotation  matrix  that  brings  I  onto  2?  is  given  by 

R1  —  iy^BIA  —  fi.g/fi.B/ )^3.v3  Bi4q BI  (166) 

Similarly,  the  desired  attitude  of  the  spacecraft  can  be  described  by  a  desired, 
body-fixed  frame  2>d  with  respect  to  7 .  In  this  manner,  the  attitude  is  parameterized  by 

considering  a  desired  unit  quaternion  c^BjI{t)  =  (t),qBdl4  (t)}  e  R3xR  where  the 

constraint  of  (105)  holds.  Therefore,  the  3x3  rotation  matrix  that  brings  I  onto  2>d  can 
be  defined  as 


,4*7  B,. 


(167) 


and  from  (108)  and  (167),  the  3x3  rotation  matrix  that  brings  2>d  onto  2>  is  therefore 


X  =  Bri  BjrT,  =(vldBA  -q^q^)43 +2q^q^  (168) 
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where  the  quaternion  error  vector  qBiB  (t)  =  jq5 /B  (t),qB/B  4  (t)|  e  M3xM  is  defined  by 
Behai  et  al.  (2002)  as 


^BdB  ^BjIA^BI  tfBjIA^BjI  +  ^BI^BjI 
QBjB  =  ^BIA^BjIA  +  ^BI^BjI 


(169) 


The  error  quaternion  kinematic  differential  equation  then  becomes 

^BdB  —  ^^BdBA^3x3  ^BjB  )  ^ BI 

'  =_I  r 

tfBjB  2  QbjB^BI 


(170) 


with  (oBI  =  (oB/B  considering  a  rest-to-rest  maneuver  where  co/(  /  =  03tl .  For  simplicity, 

the  components  of  qBB  will  be  referred  to  as  qB  B  =  {qt,q2,q3,q4}  for  the  remainder  of 

this  chapter.  The  full  rigid  spacecraft  system  of  angular  motion  equations  considering  the 
error  quaternion  kinematics  can  be  expressed  as 


Q BjB  2  (  9 8rfS, 4^3x3  +  ^BdB  )  01 BI 

■  =_I 

QBjB  2  ^BjB^BI 

til  !>[  — —J  G)BjJ(Obi  +  J  T 


(171) 


2.  Underactuated  Rigid  Spacecraft  with  Two  Control  Torques 

In  considering  an  underactuated  rigid  spacecraft  with  only  two  independent 
control  torques,  the  unactuated  axis  can  be  taken  to  be  without  loss  of  generality  and 
(96)  can  be  rewritten  as 

<*B1  =  f{™bi)™bi+t  (172) 

where  t(7)  e  M3,  Vt  =  J~‘T  =  [0, u]y  is  the  scaled  control  vector  with 
u(t)  e  U2,  Vu  =  [n2,w3]r  and  F( to)  e  IR3*3  denotes  the  drift  matrix 
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(173) 


In  this  manner,  the  full  underactuated  rigid  spacecraft  system  of  angular  motion 
equations  become 


_  1  / 

X  \ 

‘U* 

—  2  \^bjB, 4-* 3*3 

+  tfBjB  j 

1  r 

2  Q 

6)B! 

=  F((ob,)o)bi 

+  T 

(174)  can  be  written  in  control-affine  form  as 

x  =  /(x)  +  G(x)u  =  /7(x)toBy  +  Gu 


(174) 


(175) 


with  state  vector  x(t)eIR6jdR,  \/x  =  [ql,q2,q3,q4,G>l,a>2,(03]T ,  control  vector 
u(t)  e  M2,Vu  =  [m1,m2]7',  drift  matrix  F(x)elw  given  by 


F(x) 


ry  (^4^3.V3  + 


T 


—J  1  o)m  J 


(176) 


and  control  matrix  G  e  E 7-1:2  given  by 


G  = 


'5x2 


2x2 


(177) 


Furthermore,  (175)  can  be  split  into  two  coupled  differential  subsystems  of 
concern;  the  first  coupled  subsystem  combines  the  unactuated  dynamics  with  the 
unactuated  quaternion  kinematic  component,  while  the  second  considers  the  remaining 
two  quaternion  kinematic  components.  Both  of  these  systems  are  indirectly  controlled 
through  the  action  of  the  angular  velocity  components  about  the  actuated  body-fixed 


82 


axes.  The  third  subsystem  considers  the  actuated  dynamics.  In  fact,  although  global 
asymptotic  stability  cannot  be  achieved  without  unbounded  control,  the  full  system  of 
underactuated  rigid  spacecraft  equations  of  angular  motion  can  be  stabilizable  to  an 
arbitrarily  small  region  about  the  equilibrium  of  the  system  via  time-invariant  smooth 
state  feedback  control  by  considering  a  damping  coefficient  to  ensure  the  generalized 
inverse  of  the  controls  coefficient  in  the  proposed  control  law  does  not  go  unbounded  at 
times  during  the  maneuver.  This  concept  will  be  discussed  in  detail  in  subsequent 
sections. 

In  order  to  compactly  view  the  three  subsystems  of  concern,  it  is  useful  to  first 
partitioning  (173)  as  done  by  Bajodah  (2009b)  where 


F(&bi)  = 


FU  (W  BI  ) 
Fl\  ) 


FY2  (®B/  ) 
F22  (®ii  ) 


(178) 


with  Fn((dBJ)  eRlxl ,  Fl2((0BI)  eRlx2 ,F2l((0BI)  eR2xl ,  F22(o)BI)  eR2*2 ,  and  then 
consider  the  drift  matrix  associated  with  the  vector  differential  equations  of  (170)  as 


>11 1 

Fu\ 

(‘W) 

(>,«) 

F22 

(179) 


where  Fn  (q^)  e  Rlxl,  Fn  (q^)  e  RU2,F21  (q^)  e  Rlxl,  F22  (q^)  e  R2x2 .  Now,  by 

taking  the  unactuated  state  vector  to  be  xu(/)g12,Vxu  =  [qu , oju ]T  where  qu  =  qx  and 

t ou  =  cox ,  and  the  actuated  angular  velocity  vector  to  be  (oa  (7)  e  M2,  Vcafl  =  [&>2,<u3]r  ,  the 
unactuated  system  becomes 


x 


u 


0 

0 


^ii 

Fn  (®b/  ) 


Fn 


Oi 


(180) 
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By  considering  the  two  remaining  components  of  the  error  quaternion  to  be 
pseudo-actuated,  the  pseudo-actuated  error  quaternion  vector  can  be  taken  to  be 

qP(')eR2-yqP=[?2  ,q3f  and  the  actuated  error  quaternion  subsystem  then  becomes 

q „  =  [02,1  Fn  (qB(iS )]  x„  +  F22  )oa  (181) 

The  third  and  final  coupled  subsystem  which  considers  the  actuated  dynamics  and 
explicitly  contains  the  control  vector  u  is  therefore 

®a=[02xl  (182) 

B.  LINEAR  PARAMETERIZATION  OF  THE  UNDERACTUATED  RIGID 
SPACECRAFT  SUBSYSTEM 

As  mentioned  previously,  the  unactuated  subsystem  given  in  (180)  can  be  seen  to 
be  indirectly  affected  by  the  control  vector  u  through  coa .  By  defining  the  scalar 

function  0(xi():R2  ^R  as 

0(xu)=(ou+aq«  (183) 

where  a>  0  represents  a  scalar  on  qn  and  ^(xi()  is  a  continuous  twice  differentiable 
function  satisfying 

0(xu)  =  Oo6)u=-aqu  (184) 

the  unactuated  subsystem  (180)  can  now  be  transformed  into  a  stable,  linear,  second- 
order  dynamic  system  similar  to  what  is  proposed  in  Bajodah  (2009a  and  2009b)  where 

Hx« )  +  2^(x„ )  +  rV(x„ )  =  0  (1 85) 

with  the  coefficient  y  chosen  such  that  (185)  is  stable,  i.e.,  the  characteristic  equation 

s2+2ys  +  y2=0  (186) 
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has  strictly  negative-real  part  roots  thus  implying  that  y  >  0 . 


The  first  time  derivative  of  ^(xu)  along  solution  trajectories  of  the  underactuated 
rigid  spacecraft  angular  motion  equations  given  by  (175)  considering  the  control  vector 
u ,  the  control  matrix  G  given  by  (177)  andf  (x)  =  F[x)oBI  is: 

Hxu)  =  x  =  (x)  +  Gu)  =  Lf</>(xu )  +  LJ^) u  ( 1 87) 

where  LG(p{xu )  evaluates  to  the  2x1  null  vector.  Likewise,  the  second  time  derivative  of 
</>(xu)  along  the  solution  trajectories  of  the  underactuated  rigid  spacecraft  angular  motion 
equations: 

.  ,  ait^(x,).  =  3i,^(x,)(f  w  ^  (188) 

OX  OX 

where  Lfp(xu),  Ljp(xu)  andLGLfp(xu)  represent  the  Lie  derivatives,  or  directional 
derivatives  of  ^(xH)  along  the  direction  of  the  vector  fields  defined  byf(x)  and 

G(Slotine  1991).  With  (p  and  ip  given  by  (187)  and  (188),  (185)  can  then  be  expressed 
as  the  point-wise  linear  form 

ar(x)u  =  b(x)  (189) 

where  the  controls  coefficient  a(x)  e  R2xl  is  given  by 

a(x)  =  [LGLt<p(Xu)J  (190) 

and  the  scalar  controls  load  b(x)  is  given  by 

b(x)  =  -L'<p(xu)-2rLf<p(xu)-y2<p(xu)  (191) 

The  following  definitions,  propositions  and  theorem  include  only  minor  changes 

from  Bajodah  (2009b)  to  further  expand  them  to  include  the  full  state  x  containing  both 
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the  dynamics  and  kinematics  of  the  system.  The  proofs  of  these  flow  directly  from  those 
presented  in  Bajodah  (2009b)  due  to  the  cascading  nature  of  the  controls. 

Definition  1:  The  desired  linear  unactuated  dynamics  given  by  (185)  is  said  to 
be  realizable  by  the  rigid  spacecraft  angular  motion  equations  (174)  at  specific  values  of 
\(t)  if  there  exists  a  control  vector  u(t)  that  solves  (189)  for  these  values  of  \(t) .  If  this 

is  true  for  all  x(t)  such  that  xu  (t)  ^  02xl,  then  the  linear  unactuated  dynamics  given  by 
(185)  is  said  to  be  globally  realizable  by  (1 74). 

Proposition  1:  By  letting  (190)  be  the  controls  coefficient  relative  to  (183)  along 
/(x)  for  the  linear  dynamics  given  by  (185)  that  is  globally  realizable  by  (174),  then 

a(x)  =  °2u^$Kx^)  =  02u  (192) 

Proof:  For  a  vector  u  to  exist  that  solves  (189)  at  specific  values  of  x,  b(x) 
must  be  in  the  range  space  of  ar  (x)  at  that  value  of  x  .  This  is  possible  for  any  value  of 
b(x)  only  if  (189)  is  consistent.  That  is  to  say  that  not  all  elements  of  ar  (x)  vanish  at 
that  value  of  x.  Therefore,  the  existence  of  a  vector  x*  ^  0  such  that  a(x*)  =  02Tl 

implies  that  the  desired  linear  dynamics  of  (185)  are  not  realizable  at  x*  ,  which  in  turn 
violates  global  realizability  of  these  desired  linear  dynamics  and  thus  proves  sufficiency. 
Necessity  flows  from  the  fact  that  the  elements  of  the  drift  vector  f  (x)  are  multivariable 

polynomials  with  only  quadratic  elements  in  the  components  of  x  and  that  ^(xi()  has  a 

bounded  xi;  gradient,  so  that  a(06xl)  =  02xl. 

Definition  2:  The  zero  actuated  full-state  Jacobian  of  the  controls 
coefficient  (x)  is  the  square  matrix  that  results  by  partially  differentiating  the  controls 

coefficient  with  respect  to  the  actuated  angular  rates  0)a  evaluated  at  toa  =  02rl . 


/oW 


5ar (x) 


(193) 
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Proposition  2:  The  unactuated  linear  dynamics  given  by  (185)  is  globally 
realizable  by  the  rigid  spacecraft  angular  motion  equations  given  in  (1 74)  if  and  only  if 


det[yo(x)]^0 


(194) 


Proof:  The  partial  derivative  of  (189)  with  respect  to  oa  and  evaluating  the 
resulting  equation  at  o)a  =  0  yields 


y'o(x)u  = 


<96  (x) 

do.),. 


JtO  =0, 


(195) 


Invertibility  of  the  zero  actuated  state  Jacobian  of  the  controls  coefficient  implies 
that  a  control  law  ucan  be  derived  for  global  realizability  of  the  desired  linear  dynamics 
of (185)  as 


u  = 


/.■‘W 


56  (x) 


do. 


JCD  =0? 


(196) 


which  proves  necessity.  Consider  the  non-linear  time-varying  system  given  by  the 
equations 

z0=a(z0)  (197) 


where  z0  =  [zH,za]r ,  \/zu  eM,Vza  el2'1.  From  Proposition  1,  global  realizability  of  the 
desired  linear  dynamics  of  (185)  implies  that  the  origin  za=02xl  is  the  unique 
equilibrium  point  of  (197)  at  zu  =  0  .  Furthermore,  since  a(z0)  is  a  smooth  vector  field, 
it  follows  from  Milnor’s  theorem  (Bajodah  2006b)  that  it  is  also  a  global  diffeomorphism 
on  M6*1 ,  i.e.  that  it  has  continuous  partial  derivatives  and  an  invertible  zero  actuated  state 
Jacobian  and  thus  sufficiency  is  obtained. 

Theorem  1:  Given  a  non-singular  zero  actuated  state  Jacobian  J{]  ( x )  of  the 
controls  coefficient  a  (x)  =  02xl  along  /  ( x )  =  6'  ( x )  co  for  all  x(7)  eR6 ,  the  infinite  set  of 
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control  laws  that  globally  realize  the  unactuated  dynamics  by  the  underactuated  rigid 
spacecraft  angular  motion  equations  given  by  (174)  is  given  by 

u  =  u  +  P(x)y  (198) 

where  u  e  1R2  is 

u  =  a+(x)b(x)  (199) 


and  a'  (x)  e  R2aI  represents  the  Moore-Penrose  or  generalized  inverse  of  the  controls 
coefficient a(x)  e  R2'1  so  that 


(9= 


a 


(9 


K9ir 

®2xl 


|a(x)H° 

|a(x)|  =  ° 


(200) 


and  P(s.)  el 'f  represents  the  null-space  projection  matrix  of  a  (if)  such  that 


/5(x)  =  /2i2+a+(x)ar(x) 


(201) 


and  y  e  M2xl  is  an  arbitrarily  selected  null-control  vector. 

Proof:  Satisfaction  of  condition  (194)  implies  that  the  desired  linear  dynamics  of 
(185)  are  globally  realizable  by  the  underactuated  rigid  spacecraft  system  of  angular 
motion  equations.  From  Proposition  1,  this  global  realizability  further  implies  that 
a(x)  ^  02jd  at  which  infinite  number  of  solutions  for  the  point- wise  linear  relation  (189) 

exist.  Multiplying  both  sides  of  (198)  by  ar  (x)  yields 


ar  (x)u  =  ar  (x)a+  (x)b(x)  +  ar  (x)P(x)y 

=Mx)  ( 

and  thus  (189)  is  recovered.  Therefore,  the  control  vector  u  given  by  (198)  linearly 
parameterizes  all  solutions  of  (189)  by  the  null-control  vector  y .  □ 
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In  this  manner,  the  infinite  set  of  closed-loop  full  underactuated  rigid  spacecraft 
system  of  angular  motion  equations  in  partitioned  form  become 


x  = 


0  Fn(q^) 

0  Fu(o)b/)_ 


X.  + 


Fn  ( 0)  BI ) 


0) 


(203) 


®2xl  ^21  )  X«  +  ^22  )  01  a 


(204) 


94  = 


0  -\q. 


x„+qP« 


p  a 


=[°2xl  i72l(^/)]x„+F22((Oi;/)(Ofl+U 


(205) 

(206) 


(198)  consists  of  two  parts.  The  first  part,  given  by  u,  is  tenned  the  particular  solution 
and  acts  specifically  on  the  range-space  of  the  generalized  inverse  of  the  controls 
coefficienta+  (x) .  The  second  part,  given  by  P(x)y  ,  is  tenned  the  auxiliary  solution  and 
resides  in  the  orthogonal  complement  subspace,  or  the  null-space  of  the  controls 
coefficient  in  M3  with  the  null-control  vector  y  being  projected  onto  this  space  by  means 
of  the  projection  matrix  P(x) . 

As  discussed  in  Bajodah  (2009b),  the  null-control  vector  yeM2xlis  not  fully 
arbitrary  and  should  be  chosen  to  yield  stability  of  the  closed  loop  system  of  equations 
given  in  (203)  through  (206).  Two  such  null-control  vectors  will  be  presented  in  the 
following  sections,  the  first  is  Lyapunov  function  based  to  yield  global  stability  for  the 
system  given  established  bounds  on  the  control  and  the  second  is  perturbed  feedback- 
linearizing  to  yield  local  stability  and  guaranteed  boundedness.  The  development  of 
these  null-control  vectors  and  their  subsequent  proofs  of  stability  provide  a  key 
contribution  to  the  body  of  knowledge  with  respect  to  underactuated  rigid  body  attitude 
control. 
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C.  DEVELOPMENT  OF  THE  QUATERNION  FEEDBACK  REGULATOR 

AND  STABILITY  ANALYSIS 

In  this  section,  the  stability  of  the  closed-loop  system  of  equations  given  by  (203) 
through  (206)  is  addressed  for  two  derived  null-control  vectors.  From  the  definition  of 
the  generalized  inverse  of  a(x)  given  by  (200), 

lim  a+  (x)  =  co2xl  (207) 

which  implies  a  singularity  on  the  closed-loop  stability  of  the  system.  Similarly  from 
Proposition  1,  if  the  linear  unactuated  dynamics  (185)  are  globally  realizable  by  the 
underactuated  rigid  spacecraft  system  of  angular  motion  equations  given  in  (175),  then 

JT  na(X)  =  0M  (208) 

In  order  to  properly  bound  the  control  input  to  allow  for  the  control  law 
derivations  that  follow,  the  damped  generalized  inverse  and  its  corresponding  damped 
controls  coefficient  null-projector  are  used  as  presented  in  Bajodah  (2009a).  The  damped 
controls  coefficient  generalized  inverse  is  fonnulated  by  considering  an  arbitrarily  small 
damping  coefficient  /?,  to  provide  a  bound  on  the  generalized  inverse  as  the  squared 
Euclidean  nonn  of  the  controls  coefficient  tends  to  zero.  Thus,  the  damped  generalized 
inverse  of  the  controls  coefficient  a^  (x)  e  R2  vl  is 


'  a(x) 

II  /  \l|2 


l  A2 


la(x)IM 

|a(x)||<A 


(209) 


where  scalar  /)  is  a  positive  damping  coefficient.  Bounding  the  generalized  inverse  of 
a(x)  given  by  (200)  in  this  manner  smoothes  the  infinite  set  of  control  laws  presented  in 
Theorem  1.  Furthermore, 
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(210) 


and 


lim  ||a^(x)||  =  0  (211) 


and  aj(x)  pointwise  converges  to  a+(x)  as  /(vanishes  (Bajodah  2009a). 

Likewise,  the  damped  null-projection  matrix  of  a(x)is  modified  from  (209)  such 

that 


^(X)  =  /2,2-<(X)ar(X) 


(209)  and  (212)  imply  that 


pd  W = 


2x2 


l(x)|' 


||a(X)|<A 


ft 


(212) 


(213) 


and  from  (208),  during  the  steady-state  phase  of  response  the  damped  null-projection 
matrix  of  a (x)  becomes 


Jim  A(x)  =  /2,2  (214) 

such  that  the  auxiliary  part  of  the  infinite  set  of  control  laws  expressed  in  (198)  converge 
to  the  null-control  vector  y.  Construction  of  these  control  laws  with  (209)  and  (212) 

yields  the  damped  control  vector  u  defined  as 

»/  =  »rf+^(x)y  (2i5) 


where  urf  e  R2xl  is  given  by 
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(216) 


=  <(X)MX) 

Furthermore,  the  damped  control  vector  yields  the  closed-loop  full  underactuated 
rigid  spacecraft  system  of  angular  motion  equations 


x  = 


0  Fn(q^) 

0  Fn(o}5/)_ 


X  + 


Fn 

Fn  (o)/;/ ) 


Oi 


^2x1  A I  )  X«  +  ^22  ( *\-BdB  )  01 a 


9 4  = 


0 


x„+q 


[®2.vl  Al  (“5/ )]  x„  +  ^22  (®  i/ )  ®a + urf + A  (x)  y 


(217) 


(218) 

(219) 

(220) 


D.  LYAPUNOV  FUNCTION-BASED  QUATERNION  FEEDBACK 
REGULATOR 

Let  the  null-control  vector  y  be  chosen  as 


y  =  (x)[-<*>  %  (®*/ )  -  «I^2  (®« )®B  -  -  *q*,*«>]  -  (22 1) 


where  k,d  are  positive  scalar  gains  and  the  damped  actuated  dynamics 
projector  q)  (x)el2  is  defined  as 


qUx)  = 


O) 


°>Ia(x)°l 

A 


«lA(x)o>a<A 


(222) 


where  /f  is  a  positive  damping  coefficient  on  q  (x ) .  In  this  manner,  substitution  of  (221) 
into  (212)  yields  the  class  of  control  laws 

=»„  +  A(x){<(x)[-«lAi(«*/H  -<A2  (««/)«*  -<u-Aq^a>] -</©„}  (223) 
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which  yield  the  closed-loop  dynamical  actuated  subsystem 


®a  R2\  mu  R22  (WB/  )®a  +  Urf 

-0^21  (®*/  )  G>u  -  ™TaR 22  (®fl/  )  «>« 


<(x) 


T 

1  bhbP*bi 


-do)r, 


(224) 


Conjecture  1.  Let  the  function  0(xu)  be  globally  twice  continuously 
differentiable  and  satisfy  the  condition  given  by  (184),  the  controls  coefficient  a(x)  given 
by  (190)  and  controls  load  b(x)  given  by  (191)  relative  to  ^(xu)  of  the  linear  unactuated 
second-order  dynamics  given  by  (185)  along  the  drift  vector  f  (x)  =  F(x)c<)Hl  of  the 
underactuated  rigid  spacecraft  system  of  angular  motion  equations  given  by  (175). 
Furthermore,  let  the  damped  generalized  inverse  of  the  controls  coefficient  a)  (x)  be 

given  by  (209),  the  damped  null-projector  of  the  controls  coefficient  Pd(x)  be  given  by 
(212),  and  the  null-control  vector  y  be  given  by  (221)  with  damped  actuated  dynamics 
projector  vfd  ( x )  be  given  by  (222).  If  the  zero-actuated  state  Jacobian  ( x )  given  by 
(193)  is  non-singular,  then  the  positive  scalar  feedback  gains  k,d  e  M  can  be  chosen  to 
globally  stabilize  the  closed-loop  underactuated  rigid  spacecraft  system  of  angular 
motion  equations  given  by  (217)  through  (220)  by  bounded,  smooth,  time-invariant  state 
feedback  to  an  arbitrarily  small  vicinity  of  the  equilibrium  of  the  system  defined  by  the 
damping  coefficients  f:(  and  fi2 . 

Proof.  Let  the  scalar  function  e(x)  be  defined  as 

e(x)  =  ^(x„)  +  ^(x„)  (225) 

where  y  is  as  defined  for  (186).  Given  (189)  and  (225)  ,  e(x)  therefore  evaluates  to 

e(x)  =  ^(xH)  +  ^(x„)  =  -ye(x)-Z)(x)  +  ar(x)u  (226) 

Let  the  Lyapunov  control  function  be  defined  as 
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(227) 


v  =  \e2  (x)  +  V>„  +  7, 2  +  +  («4  -  lf 

=  ^2(x)+^r‘toX+2(i-^4) 

Note  that  V  is  positive  definite  and  asymptotically  unbounded  in  e(\)  andoy , . 

Furthermore,  V  can  be  partitioned  into  two  positive  definite  functions  such  that 
V  =  Vl  +  V2  where 


^l  =  ^e2(x) 

(228) 

^X+2(l-?4) 

(229) 

where  k  1  e  R  is  a  positive  scalar  constant  gain. 

The  corresponding  time-derivative  of  Vx  is  given  by 

Vx  =e(x)e(x)  =  -e(x)h(x)-7e2(x)  +  e(x)ar(x)u  (230) 


which  by  is  reduced  by  substitution  of  (226)  given  that  ar  (x)Pd  (x)  =  0  by  definition  to 

V]=-re2(x)  (231) 

Note  that  (231)  is  guaranteed  to  be  negative  definite  independent  of  the  choice  on 
the  null-control  vector  y  with  proper  selection  of  y . 

Given  that  the  null-control  vector  y  can  be  arbitrarily  chosen,  it  can  be  selected  to 
ensure  V2  <  0  which,  when  combined  with  (231),  ensures  that  V  is  negative  definite  and 
global  asymptotic  stability  of  the  underactuated  rigid  spacecraft  system  is  guaranteed. 
The  corresponding  time-derivative  of  V2  is 


V2  h 

=  ^(olF2{  (<o) 0)u  +  k~l(olF22  (e»)e»a  +  k~l(oTa u  +  k~WaPd  (x) y  +  q> 


(232) 
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Let  the  null-control  vector  y  be  defined  as 


with  a  positive  scalar  gain  del.  Similar  to  the  damping  of  the  control  coefficient 
generalized  inverse,  the  damped  actuated  dynamics  generalized  inverse  vfd  (x)  given  by 
(222)  implies  that 

h;wlls4-  <234) 

Pi 

and  that 

lim  L+(x)|  =  0  (235) 

oa^02^i  H  v  7 II 

and  that  (x)pointwise  converges  to  x\+d  (x)as  /?,  vanishes. 

In  the  condition  where  to TaPd  (x)  =  0lx2,  the  null-control  vector  can  be  seen  to  not 

affect  the  stability  of  the  system  in  a  Lyapunov  sense  from  (232).  However,  do  to  the 
arbitrarily  small  nature  of  this  region  by  selection  of  the  damping  coefficient  /f ,  it  is 
conjectured  that  the  system  is  globally  stable  to  this  arbitrarily  small  region  defined  by 
the  damping  coefficients  /)  and  /L  with  the  feedback  control  law 

u  =  a^(x)h(x)  +  Prf(x)y  (236) 

where  y  is  given  by  (233). 
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E.  FEEDBACK-LINEARIZING  QUATERNION  FEEDBACK  REGULATOR 

Let  the  null-control  vector  y  be  chosen  as 

y  =  -Fn  (to)  o)u  -  F22  (to)  cofl  -  d(na  -  kqp  (237) 

where  d,k  are  positive  scalar  gains  to  be  detennined.  In  this  manner,  substitution  of 
(237)  into  (212)  yields  the  class  of  control  laws 

=Urf  +^(x)[-^21(o>H  -F22((o)o)a -d(0a  -kqp~]  (238) 

The  close-loop  actuated  dynamical  subsystem  becomes 

tba  =  F21  (to)  o)u  +  F22  (to)  tofl  +  ud  +P,(x) [~F21  (to) -  F22  (to) to„  -  d(oa  -  kqp ]  (239) 

which  converges  by  (210)  and  (214)  to  the  linear  dynamical  system  (Bajodah  2009a) 

™a=-d(»a-kqp  (240) 

Furthermore,  the  linear  dynamical  system  given  by  (240)  is  upper  bounded  by 
selection  of  the  damping  coefficient  /)  forming  a  domain  of  attraction  in  x  whereby 

||a(x)||>  f)  due  to  the  damped  controls  coefficient  generalized  inverse  defined  by  (209). 

Therefore,  the  dynamical  system  given  by  (240)  coupled  with  the  remaining  equations  of 
motion  given  by  (217)  through  (219)  form  a  feedback  linearizing  transfonnation  from  the 
underactuated  rigid  spacecraft  angular  motion  equations  given  by  (175)  over  a  domain  of 
attraction  in  x . 

Conjecture  2.  Let  the  function  ^(xu)  be  globally  twice  continuously 
differentiable  and  satisfy  the  condition  given  by  (184),  the  controls  coefficient  a(x  )  given 
by  (190)  and  controls  load  b(\)  given  by  (191)  relative  to  of  the  linear  unactuated 

second-order  dynamics  given  by  (185)  along  the  drift  vector  /(x)  =  F(\)oiljl  of  the 
underactuated  rigid  spacecraft  system  of  angular  motion  equations  given  by  (175) 
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Furthermore,  let  the  damped  generalized  inverse  of  the  controls  coefficient  (x)  be 
given  by  (209),  the  damped  null-projector  of  the  controls  coefficient  Pd(x)be  given  by 
(212),  and  the  null-control  vector  y  be  given  by  (221).  If  the  zero-actuated  state 
Jacobian  (x)  given  by  (193)  is  nonsingular,  then  the  positive  scalar  feedback  gains 
k,d  eR  can  be  chosen  to  stabilize  the  closed-loop  underactuated  rigid  spacecraft  system 
of  angular  motion  equations  given  by  (217)  through  (220)  by  bounded,  smooth,  time- 
invariant  state  feedback  control  to  an  arbitrarily  small  vicinity  of  the  equilibrium  of  the 
system  defined  by  the  damping  coefficients  f:(  within  a  domain  of  attraction  in  x  (r ) . 

Proof:  Assuming  the  inverse  damping  gain  k~x  e  R  exists  and  is  positive 
definite,  consider  the  control  Lyapunov  function 


1  2 
V  =  V,  +  +  <?2  +  9.1  +  (<?4  _  1) 

=  ^0)Tak-l(Oa+2(\-q4) 


(241) 


Note  that  V  is  positive  definite.  This  particular  Lyapunov  function  is  well  known 
for  spacecraft  attitude  dynamics  appearing  in  Wie  and  Barba  (1985);  Wie,  Weiss  and 
Arapostathis  (1989);  and  Cristi  and  Burl  (1993)  among  others. 

Differentiating  V  along  the  trajectories  of  the  underactuated  system  given  by 
(175)  yields 


V  =  (i)Takl6)a-2q4 
=  -dfk  'doa  +qu0Ju 


(241)  can  be  rewritten  as 


F  =  L||ra„f  +2(l-?4) 


(242) 


(243) 


and  (242)  can  be  rewritten  as 
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+  q«®« 


(244) 


V  =  -M||coa|" 

Given  ^(xt()  defined  by  (183)  such  that 

<!>{\)  =  0Ju+aclu 


(245) 


(244)  can  be  rewritten  as 


V  =  qu</>{*u)-aq2u-kd 


to , 


(246) 


From  Lemma  1  in  Appendix  B,  |coa  (7)||  is  uniformly  bounded, 
i.e,.  |coa  (7)||  <c\/t .  By  defining  /(t)eK  to  be  the  scalar  function  described  by 


/W=jlK(r)l" dT  (247) 

0 


and  given  the  designed  dynamics  in  (185),  ^(xu)(7)  — >  0  exponentially,  application  of 
Lemma  2  in  Appendix  B  yields 

lim/(7)  =  /(oo)<+oo  (248) 

?->+oo  V  7  V  7 

and 


f(t)  =  i||coa  (r)|2  d,T  <  c  \/t  (249) 

0 

Furthermore,  with  (247)  and  (249)  defined,  Barbalat’ s  Lemma  (Slotine  and  Li 
1991)  can  be  applied  by  conjecture  given  an  arbitrarily  small  selection  of  the  damping 
coefficient  f3x  on  the  controls  coefficient  generalized  inverse  to  yield 


lim  / (V)  =  lim  (0<;  (V)  =  0 


(250) 


From  Lemma  3  in  Appendix  B, 
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(251) 


HMOH 

and  given  Lemma  4  in  Appendix  B,  (250)  and  (25 1)  yield 

Hkwil=o  (252> 

t-¥ oo  N  y  v  7  II 

Thus,  the  closed-loop  equations  of  motion  given  by  (217)  through  (220)  are  stable 
within  a  domain  of  attraction  in  x(t) . 

F.  DESIGN  EXAMPLE 

In  order  to  demonstrate  the  effectiveness  of  the  presented  quaternion  feedback 
regulators  with  both  null-control  vector  selections  for  underactuated  spacecraft,  a  control 
design  example  is  presented  that  considers  an  asymmetric  rigid  spacecraft  with  principal 
moments  of  inertia  (in  kgm 2)  of  Jn  =32.5,  J22  =25  and  J33  =  12.5  .  These  selections  of 

moments  of  inertia  correspond  to  a  30  kg  rigid  rectangular  spacecraft  with  uniform  mass 
distribution  and  sides  of  length  1x2x3  m.  For  purposes  of  representing  real-world 
spacecraft  torqueing  capabilities,  the  maximum  available  torque  about  each  control  axis 
is  limited  to  1  Nm. 

From  (190),  the  control  coefficient  a(x)  relative  to  the  feedback  linearizing 
transformation  0(xw )  given  by  (183)  is 


a 


(0 


1 

2^ 


(J22  ~  T33 ) 
J]: 


CO , 


The  determinant  of  the  associated  (x)  as  defined  by  (193)  is 


(253) 
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=  det 


0 


(254) 


det 


daT (x) 


5oo. 


CO  =0^ 


{^22  "^33) 

Jn 


(J22  J22) 

Ju 


to  =0-, 


which  is  non-zero  for  all  x(t)  e  R7rl . 

For  the  numerical  simulations  that  follow,  a  fourth-order  Runge-Kutta  numerical 
integration  scheme  is  used  to  integrate  (217)  through  (220)  with  a  fixed  time  step  of  .Is. 
Appendix  C  provides  the  code  for  the  algorithm  implemented  via  an  Embedded 
MATLAB®  function.  In  order  to  accommodate  exact  zero  initial  conditions  on  to  and 
q  without  failure  of  the  regulator,  the  measured  state  is  perturbed  to  1x1 0  4  in  these 
situations. 

1,  Numerical  Simulations  using  Lyapunov  Function-Based  Quaternion 
Feedback  Regulator 

a.  Rest  to  Rest  Reorientation  from  Large  Initial  Angle  Offset  Given 
Low  Available  Control  Torque 

A  large-angle  rest-to-rest  reorientation  maneuver  is  considered.  The 
quaternion  describing  the  initial  orientation  (t0)of  2>d  with  respect  to  I  is  taken  to  be 

q(t0)  =  [.57,.57,.57,.159]/  which,  by  considering  a  3-2-1  sequence  of  rotations, 

corresponds  to  initial  Euler  angles  of  [cq  (t0),a2  (t0),a3  (t0)]  =  [l09.8,-27.9,109.8]deg  . 

The  final  orientation  is  taken  to  be  q^)  =  [0,0,0, l]r .  The  initial  and  final  angular 

velocity  vector  is  co(70)  =  w(ty  )  =  [0,0, 0]r .  This  reorientation  maneuver  is  chosen  in 

order  to  compare  against  the  fully  actuated  quaternion  feedback  controller  presented  in 
Wie  and  Barba  (1985).  The  chosen  gains  are 

/  =  .70,  a  =  1 .25,  d  =  7.5,  k  =  2.25,  f,f2  =  lxl  O'4 
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Figure  15  shows  the  time  histories  of  the  Euler  angles  and  angular  rates  and 
demonstrates  smooth  attitude  stabilization  within  90  s.  Figure  16  shows  the  time 
histories  of  the  control  torques.  The  control  noise  near  the  inherent  singularity  in  the 
null-control  vector  given  by  (233)  is  evident  throughout  the  maneuver  until  the 
stabilization  region  is  reached. 
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TIME  (SEC) 


Figure  15.  Time  Histories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Rest  to 
Rest  Reorientation  Manuever  from  Large  Initial  Angle  Offset  Using  the  Lyapunov- 
Function  Based  Control  Law  and  Nominal  Available  Control  Torque 
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Figure  16.  Time  Histories  of  the  Control  Torques  for  a  Rest  to  Rest  Reorientation 
Manuever  from  Large  Initial  Angle  Offset  Using  the  Lyapunov-Function  Based  Control 
Law  and  Nominal  Available  Control  Torque 
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b.  Large  Initial  Angular  Rate  to  Rest  Stabilization  and 
Reorientation  Given  Low  Available  Control  Torque 

A  three-axis  detumbling  and  reorientation  maneuver  is  considered.  The 
quaternions  describing  the  initial  and  final  orientation  of  2>d  with  respect  to  I  are 

=  q(f0)  =  [0,0,0, l]r and  the  initial  angular  velocity  vector  is  taken  to  be 

co(t0)  =  [l,-l,l]ra<af /s2  .  The  final  angular  velocity  vector  is  «(?/)  =  [0,0,0]r  .  By 

assuming  the  same  choice  of  parameters  given  above,  Figure  17  depicts  the  time  histories 
of  the  Euler  angles  and  angular  rates  and  demonstrates  attitude  stabilization  within  300  s. 
The  failure  to  smoothly  converge  to  the  origin  is  due  to  the  limited  control  torque.  Figure 
18  shows  the  time  histories  of  the  control  torques  and  again  shows  significant  control 
effort  noise  during  the  maneuver  until  the  stabilization  region  is  reached. 
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Figure  17.  Time  Histories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Large 
Initial  Angular  Rate  to  Rest  Reorientation  Manuever  Using  the  Lyapunov-Function 
Based  Control  Law  and  Nominal  Available  Control  Torque 
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Figure  18.  Time  Histories  of  the  Control  Torques  for  a  Large  Initial  Angular  Rate  to 
Rest  Reorientation  Manuever  Using  the  Lyapunov-Function  Based  Control  Law  and 

Nominal  Available  Control  Torque 
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c.  Large  Initial  Angular  Rate  to  Rest  Stabilization  and 
Reorientation  Given  Large  Available  Control  Torque 

A  three-axis  detumbling  and  reorientation  maneuver  is  considered  again 
with  increased  torque  capabilities  as  is  typically  studied  in  most  literature  involving 
underactuated  spacecraft  control.  The  quaternions  describing  the  initial  and  final 

orientation  of  2?d  with  respect  to  I  arc  again  taken  to  be  q  [tf  j  =  q  (t0 )  =  [0, 0, 0,  l]  and 
the  initial  angular  velocity  vector  is  again  taken  to  be  co(t0)  =  [l,-l,l]rac/  /  s2  .  The  final 
angular  velocity  vector  is  to [tf )  =  [0, 0, 0]r  .  The  maximum  torque  is  taken  to  equal  to  the 

moment  of  inertia  about  the  second  principal  axis  so  that  u2  is  limited  to  1  rad  Is 2 .  By 

assuming  the  same  choice  of  parameters  given  above,  Figure  19  shows  the  time  histories 
of  the  Euler  angles  and  angular  rates,  demonstrating  attitude  stabilization  in  less  than  90s. 
Figure  20  shows  the  time  histories  of  the  control  torques.  The  control  noise  is  clearly 
evident  near  the  inherent  singularities  in  the  null-control  vector  and  its  chattering  effect 
on  the  angular  rates. 
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Figure  19.  Time  Histories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Large 
Initial  Angular  Rate  to  Rest  Reorientation  Manuever  Using  the  Lyapunov-Function 
Based  Control  Law  and  Large  Available  Control  Torque 
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Figure  20.  Time  Histories  of  the  Control  Torques  for  a  Large  Initial  Angular  Rate  to 
Rest  Reorientation  Manuever  Using  the  Lyapunov-Function  Based  Control  Law  and 

Large  Available  Control  Torque 
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d.  Rest  to  Rest  Reorientation  Manuever  Given  an  Initial  Angle 
Offset  about  Only  the  Unactuated  Axis  and  Low  Available 
Control  Torque 

An  attitude  error  stabilization  maneuver  is  considered  about  the 
unactuated  body-fixed  axis  with  low  torque  capabilities.  The  quaternion  describing  the 

initial  orientation  (t0)of  21  d  with  respect  to  I  is  taken  to  be  q(t0)  =  [,17,0,0,.98]r  which, 
by  considering  a  3-2-1  sequence  of  rotations,  corresponds  to  initial  Euler  angles  of 
[cq  (t0),a2  (t0),a3  (t0)]  =  [20, 0,0] deg .  The  final  orientation  is  taken  to  be 

q(^)  =  [0,0, 0,1]  .  The  initial  and  final  angular  velocity  vector  is 

co(t0)  =  to(^)  =  [0, 0, 0]r  .  By  assuming  the  same  choice  of  parameters  given  above, 

Figure  2 1  shows  the  time  histories  of  the  Euler  angles  and  angular  rates,  demonstrating 
attitude  stabilization  in  less  than  70  s.  Figure  22  shows  the  time  histories  of  the  control 
torques  and  demonstrates  significant  control  noise  near  the  singularities  in  the  control 
algorithm  despite  the  small  initial  angular  offset. 
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Figure  2 1 .  Time  Flistories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Rest  to 
Rest  Reorientation  Manuever  Given  an  Initial  Angle  Offset  about  Only  the  Unactuated 
Axis  Using  the  Lyapunov-Function  Based  Control  Law  and  Nominal  Available  Control 

Torque 

111 


1 


0.8 


0.6 


0.4 


0.2 

E 

5-  0 

CN 

-0.2 


-0.4 


-0.6 


-0.8 


-1 


0  20  40  60  80  100  120  140  160  180 

TIME  (SEC) 


1 


0.8 


0.6 


0.4 


0.2 

E 

5-  0 
CO 

-0.2 


-0.4 


-0.6 


-0.8 


-1 


0  20  40  60  80  100  120  140  160  180 

TIME  (SEC) 

Figure  22.  Time  Histories  of  Control  Torques  for  a  Rest-To-Rest  Reorientation 
Manuever  Given  an  Initial  Angle  Offset  about  Only  the  Unactuated  Axis  Using  the 
Lyapunov-Function  Based  Control  Law  and  Nominal  Available  Control  Torque 
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e.  Attitude  Maintenance  in  the  Presence  of  Relatively  Large 
Disturbance  Torques 

In  order  to  demonstrate  a  limitation  of  the  presented  Lyapunov-function 
based  quaternion  feedback  regulator  with  respect  to  attitude  maintenance,  an  attitude 
error  maintenance  scenario  is  considered  with  low  torque  capabilities  in  the  presence  of 
relatively  large  disturbance  torques.  The  quaternion  describing  the  initial  orientation 

(t0 )  of  2>d  with  respect  to  I  is  taken  to  be  q  (/„ )  =  [0, 0, 0,  l]r  while  the  final  orientation  is 
taken  to  be  q(/y)  =  [0,0,0,  l]  .  The  initial  and  final  angular  velocity  vectors  are 

cd(70)  =  =  [0, 0, 0]r  .  The  same  choices  of  parameters  given  above  are  used  to 

include  a  maximum  torque  about  the  two  control  axes  of  1  Nm.  Furthennore,  a  relatively 
large  disturbance  torque  of  .1  Nm  is  imparted  about  each  body  axis.  For  the  presented 
rigid  spacecraft  inertia  matrix  of  J  =  diag (32.5, 25, 12.5) kgm2 ,  this  corresponds  to 

angular  accelerations  of  .003  m/ s2 ,  .004  m/ s2  and  .008  m/ s2  about  the  respective  body 
axes.  Figure  23  shows  the  Euler  angles  and  angular  rates  and  demonstrates  that  although 
the  desired  attitude  is  unable  to  be  maintained,  it  does  remain  bounded.  Figure  24  reports 
the  control  history  for  the  attitude  maintenance  scenario  and  again  shows  the  large 
control  effort  noise  given  the  near-singularities  that  exist  in  the  controller. 
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Figure  23.  Time  Histories  of  Euler  Angles  and  Angular  Velocities  for  an  Attitude 
Maintenance  Scenario  in  the  Presence  of  Large  Disturbance  Torques  Using  the 
Lyapunov-Function  Based  Control  Law  and  Nominal  Available  Control  Torque 
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Figure  24.  Time  Histories  of  the  Control  Torques  for  an  Attitude  Maintenance 
Scenario  in  the  Presence  of  Large  Disturbance  Torques  Using  the  Lyapunov-Function 
Based  Control  Law  and  Nominal  Available  Control  Torque 
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2.  Numerical  Simulations  using  Feedback  Linearized  Based  Quaternion 

Feedback  Regulator 

a.  Rest  to  Rest  Reorientation  from  Large  Initial  Angle  Offset  and 
Low  Available  Control  Torque 

A  large-angle  rest-to-rest  reorientation  maneuver  is  considered.  The  initial 
and  final  conditions  of  the  spacecraft  are  taken  to  be  equivalent  to  the  considered 
maneuver  above  with  the  Lyapunov-based  approach.  The  chosen  gains  are 

Y  =  HO, a  =  1.25, d  =  7.5, k  =  2.25,  f  =  WO  4 

Figure  25  shows  the  time  histories  of  the  Euler  angles  and  angular  rates, 
demonstrating  attitude  stabilization  within  100  s.  The  discontinuity  in  the  Euler  angles  is 
due  to  singularities  in  the  Euler  rotation  sequence  and  does  not  affect  the  algorithm. 
Figure  26  shows  the  time  histories  of  the  control  torques  and  demonstrates  a  significant 
advantage  over  the  Lyapunov-function  based  controller  with  no  control  noise. 
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Figure  25.  Time  Histories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Rest  to 
Rest  Reorientation  Manuever  from  Large  Initial  Angle  Offset  Using  the  Feedback 
Linearizing  Control  Law  and  Nominal  Available  Control  Torque 
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Figure  26.  Time  Histories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Rest  to 
Rest  Reorientation  Manuever  from  Large  Initial  Angle  Offset  Using  the  Feedback 
Linearizing  Control  Law  and  Nominal  Available  Control  Torque 
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b. 


Rate  Stabilization  and  Reorientation  (Low  Torque) 


A  three-axis  detumbling  and  reorientation  maneuver  is  considered.  The 
initial  and  final  conditions  of  the  spacecraft  are  taken  to  be  equivalent  to  the  considered 
maneuver  above  with  the  Lyapunov-based  approach.  By  assuming  the  same  choice  of 
parameters,  Figure  27  depicts  the  time  histories  of  the  Euler  angles  and  angular  rates, 
demonstrating  convergence  to  the  zero  state  in  less  than  250  s.  Figure  28  shows  the  time 
histories  of  the  control  torques  and  again  shows  significant  control  effort  savings  over  the 
Lyapunov-based  approach  without  sacrificing  realizability  for  the  system. 
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Figure  27.  Time  Histories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Large 
Initial  Angular  Rate  to  Rest  Reorientation  Manuever  Using  the  Feedback  Linearing 
Control  Law  and  Nominal  Available  Control  Torque 
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Figure  28.  Time  Histories  of  the  Control  Torques  for  a  Large  Initial  Angular  Rate  to 
Rest  Reorientation  Manuever  Using  the  Feedback  Linearizing  Control  Law  and  Nominal 

Available  Control  Torque 
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c.  Rate  Stabilization  and  Reorientation  (Large  Torque) 

A  three-axis  detumbling  and  reorientation  maneuver  is  considered  again 
with  increased  torque  capabilities  as  is  discussed  for  the  Lyapunov-based  approach.  The 
initial  and  final  conditions  of  the  spacecraft  are  taken  to  be  equivalent  to  the  considered 
maneuver  above  with  the  Lyapunov-based  approach.  By  assuming  the  same  choice  of 
parameters  given  above  to  include  a  maximum  torque  about  the  second  and  third  body 
axes  so  that  u-,  is  limited  to  1  rad  /  s2 ,  Figure  29  shows  the  time  histories  of  the  Euler 

angles  and  angular  rates.  Due  to  the  large  available  torque,  the  spacecraft  is  stabilized  to 
the  desired  attitude  in  less  than  100  s.  Figure  30  shows  the  time  histories  of  the  control 
torques  and  shows  that  even  with  a  large  available  torque,  the  control  does  not  saturate 
during  the  maneuver  as  is  the  case  with  the  Lyapunov-based  solution. 
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Figure  29.  Time  Histories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Large 
Initial  Angular  Rate  to  Rest  Reorientation  Manuever  Using  the  Feedback  Linearizing 
Control  Law  and  Large  Available  Control  Torque 


123 


25 


20 

15 


-20 

-25 


0  20  40  60  80  100  120  140  160  180 

TIME  (SEC) 


TIME  (SEC) 

Figure  30.  Time  Histories  of  the  Control  Torques  for  a  Large  Initial  Angular  Rate  to 
Rest  Reorientation  Manuever  Using  the  Feedback  Linearizing  Control  Law  and  Large 

Available  Control  Torque 
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d.  Unactuated  Attitude  Error  Stabilization  (Low  Torque) 

An  attitude  error  stabilization  maneuver  is  considered  about  the 


unactuated  body-fixed  axis  with  low  torque  capabilities.  The  initial  and  final  conditions 
of  the  spacecraft  are  taken  to  be  equivalent  to  the  considered  maneuver  above  with  the 
Lyapunov-based  approach.  By  assuming  the  same  choice  of  parameters  given  above, 
Figure  31  shows  the  time  histories  of  the  Euler  angles  and  angular  rates,  demonstrating 
smooth  reorientation  to  the  desired  attitude  in  less  than  80  s.  Figure  32  shows  the  time 
histories  of  the  control  torques  and  shows  the  advantage  in  terms  of  control  effort  over 
the  Lyapunov-based  solution. 
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Figure  3 1 .  Time  Flistories  of  the  Euler  Angles  and  Angular  Velocities  for  a  Rest  to 
Rest  Reorientation  Manuever  Given  an  Initial  Angle  Offset  About  Only  the  Unactuated 
Axis  Using  the  Feedback  Linearizing  Control  Law  and  Nominal  Available  Control 

Torque 
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Figure  32.  Time  Histories  of  the  Control  Torques  for  a  Rest  to  Rest  Reorientation 
Manuever  Given  an  Initial  Angle  Offset  About  Only  the  Unactuated  Axis  Using  the 
Feedback  Linearizing  Control  Law  and  Nominal  Available  Control  Torque 
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e.  Attitude  Maintenance  in  the  Presence  of  Relatively  Large 
Disturbance  Torques 

An  attitude  error  maintenance  scenario  is  considered  with  low  torque 
capabilities  in  the  presence  of  relatively  large  disturbance  torques  using  the  presented 
feedback  linearizing  control  law.  By  assuming  the  same  initial  and  final  conditions  and 
parameters  used  for  the  Lyapunov-based  approach,  Figure  33  shows  the  Euler  angles  and 
angular  rates  and  demonstrates  that  although  the  desired  attitude  is  unable  to  be 
maintained,  it  remains  bounded.  By  comparison  to  Figure  23  for  the  Lyapunov-function 
based  control  law,  the  feedback  linearizing  control  law  provides  an  additional  measure  of 
stability  by  affecting  a  strictly  periodic  nature  to  the  spacecraft’s  attitude  and  a  tighter 
bounds  on  the  attitude  error  departures.  Figure  34  reports  the  control  history  for  the 
attitude  maintenance  scenario  and  depicts  a  relatively  smooth  control  effort  throughout 
the  maneuver. 
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Figure  33.  Time  Histories  of  Euler  Angles  and  Angular  Velocities  for  an  Attitude 
Maintenance  Scenario  in  the  Presence  of  Large  Disturbance  Torques  Using  the  Feedback 
Linearizing  Control  Law  and  Nominal  Available  Control  Torque 
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Figure  34.  Time  Histories  of  the  Control  Torques  for  an  Attitude  Maintenance 
Scenario  in  the  Presence  of  Large  Disturbance  Torques  Using  the  Feedback  Linearizing 
Control  Law  and  Nominal  Available  Control  Torque 
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y. 


CONCLUSION 


A.  SUMMARY  OF  CONTRIBUTIONS 

A  planar  laboratory  testbed  is  first  introduced  for  the  simulation  of  autonomous 
proximity  maneuvers  of  uniquely  control  actuator  configured  spacecraft.  This  testbed 
consists  of  a  floating  robotic  simulator  via  planar  air  bearings  on  a  flat  floor  that  is 
equipped  with  dual  vectorable  cold-gas  thrusters  and  miniature  control  moment 
gyroscope.  Inertial  position  and  attitude  measurements  are  obtained  and  fused  with  a 
discrete  Kalman  filter  and  linear  quadratic  estimator  for  navigation  while  feedback 
linearized  control  coupled  with  a  linear  quadratic  regulator  and  Schmitt  trigger  logic 
directly  command  the  control  moment  gyro  and  command  dual  in-plane  vectorable  cold- 
gas  thrusters  through  a  Pulse  Width  Modulation  scheme. 

The  presented  experimental  tests  of  autonomous  closed  path  proximity  maneuvers 
of  the  spacecraft  simulator  offer  significant  sample  cases.  The  experimental  results, 
which  show  good  repeatability  and  robustness  against  disturbance  and  sensor  noise, 
validate  the  proposed  estimation  and  control  approaches  and  validate,  in  particular,  the 
analytical  small-time  local  controllability  of  the  uniquely  configured  system.  The 
achieved  accuracy  in  following  the  reference  trajectory  (respectively,  ~  1  cm  for 
translation  and  ~  .5  deg  for  rotation  given  only  the  vectorable  thrusters  as  control  inputs) 
demonstrates  both  a  feasible  and  promising  actuator  configuration  for  small  spacecraft 
design  given  size,  weight  and  fuel  storage  considerations. 

Leveraging  the  research  on  the  3 -DoF  robotic  testbed,  a  rigorous  small-time  local 
controllability  analysis  is  conducted  to  detennine  the  minimum  number  of  control 
actuators  for  full-order  spacecraft  in  a  proximity  operations  environment.  Two  potential 
configurations  flow  from  this  analysis.  The  first  involves  a  traditional  approach  of  fixed 
thrusters  whereby  two  fixed  thrusters  are  oppositely  mounted  and  used  to  translate  the 
spacecraft  while  two  sets  of  paired  thrusters  provide  attitude  control.  The  second 
configuration  involves  a  novel  design  consisting  of  only  two  hemispherically  vectorable 
thrusters.  By  combining  the  vectorable  thrusters  in  an  opposing  manner,  the  typical 
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spacecraft  system  requiring  six  or  more  control  actuators  to  achieve  6-DoF  control  is 
reduced  to  only  two.  Furthermore,  by  simply  adding  two  additional  paired  thrusters  to 
provide  a  torque  about  the  vector  that  joins  the  two  vectorable  thrusters,  this  system 
becomes  feedback  linearizable  and  similar  methods  to  those  presented  for  the  3 -DoF 
robotic  testbed  can  be  employed  to  control  both  its  attitude  and  position. 

By  specifically  focusing  on  the  minimal  number  of  control  actuators  for  the  6- 
DoF  system,  that  is  to  say  where  only  two  vectorable  thrusters  are  required  to  provide 
small-time  local  controllability,  a  control  problem  exists  where  one  must  achieve  three- 
axis  stabilization  of  the  attitude  with  only  two  control  torques.  To  solve  this  problem,  a 
novel  quaternion  feedback  regulator  is  presented  that  capitalizes  on  recent  developments 
in  generalized  inversion  and  perturbed  feedback  linearizing  control.  A  desired  second- 
order  linear  dynamics  in  a  function  of  the  angular  velocity  and  error  quaternion 
components  about  the  unactuated  axis  is  evaluated  along  trajectories  of  the  3 -DoF 
spacecraft  angular  motion  equations  that  include  both  Euler’s  dynamic  equations  and  the 
error  quaternion  kinematics.  The  control  variables,  which  are  composed  of  a  particular 
and  auxiliary  part,  are  detennined  to  yield  a  stable  underactuated  system.  The  particular 
part  seeks  to  realize  the  desired  linear  dynamics  through  generalized  inversion  of  the 
controls  coefficient  while  the  null-control  vector  in  the  auxiliary  part  is  selected  to  yield 
stability  of  the  full  underactuated  system.  Two  control  design  methodologies  are  derived 
that  involve  separate  constructions  of  the  null-control  vector.  The  first  is  Lyapunov 
based  and  yields  a  stable  underactuated  spacecraft  system  while  the  second  involves 
perturbed  feedback  linearization  and  yields  stability  of  the  system  within  a  domain  of 
attraction.  For  both  cases,  in  order  to  overcome  the  potential  for  the  control  laws  to 
produce  numerical  instability  as  the  generalized  inverse  of  the  controls  coefficient 
becomes  singular,  the  generalized  inverse  is  damped  near  a  bound  on  the  singularity,  thus 
providing  smooth  control  and  resulting  in  stabilization  near  to  the  zero  error  state.  As 
such,  neither  controller  contradicts  the  conjecture  by  Byrnes  and  Isidori  (1991)  that  no 
smooth,  time-invariant  state  feedback  controller  can  be  found  to  locally  asymptotically 
drive  the  system  to  the  zero  error  state.  However,  given  the  arbitrary  selection  of 
damping  coefficients,  the  region  about  the  zero  error  state  can  be  reduced  to  an 
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insignificant  sized  ball.  Several  simulations  results  are  presented  for  an  underactuated 
asymmetric  spacecraft  with  two,  realistically  bounded,  body-fixed  torques  to  demonstrate 
the  wide-reaching  capabilities  of  both  presented  null-control  vector  designs.  Included  are 
a  large-angle  rest-to-rest  reorientation  maneuver,  a  large  angular  velocity  detumbling  and 
simultaneous  reorientation  maneuver  and  a  single-axis  reorientation  about  the  unactuated 
axis.  The  proposed  attitude  control  method  is  not  intended  to  provide  attitude 
maintenance  or  for  attitude  tracking  in  the  presence  of  relatively  large  disturbance 
torques;  however,  it  may  prove  widely  applicable  to  detumbling  and  reorientation 
maneuvers  of  spacecraft  with  only  two  available  control  torques. 

B.  POSSIBLE  FUTURE  DEVELOPMENTS 

The  following  iterating  research  remains  open: 

1.  Validate  by  numerical  and/or  experimental  methods  the  feedback  linearizing 
control  method  proposed  in  Chapter  III  for  the  full  6-DoF  relative  dynamics 
problem. 

2.  Determine  a  roto-translation  control  strategy  for  a  spacecraft  with  two 
oppositely  mounted  hemispherical  thrusters.  A  possible  solution  is  to 
investigate  using  a  Frenet  frame  waypoint  tracking  solution  where  the  forward 
face  is  first  made  normal  to  the  path  by  one  of  the  quaternion  feedback 
regulators  presented  in  Chapter  IV  and  then  a  separate  controller  is  used  to 
translate  the  spacecraft. 

3.  Provide  a  formal  proof  of  global  stability  for  both  the  Lyapunov-function 
based  and  feedback  linearizing  underactuated  quaternion  feedback  controllers 
presented  in  Chapter  IV.  Given  the  presented  numerical  simulation  results 
that  demonstrate  a  very  large  domain  of  attraction,  it  is  possible  that  global 
stability  for  both  controllers  can  be  rigorously  demonstrated. 

4.  Compare  both  underactuated  attitude  controllers  to  the  optimal  solution  for 
each  presented  case.  In  this  manner,  a  method  of  gain  tuning  may  be  achieved 
to  ensure  a  near-optimal  solution  can  be  achieved. 
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APPENDIX  A.  MATLAB  CODE  TO  PERFORM  SMALL  TIME 
LOCAL  CONTROLLABILITY  STUDIES 


function  [L,  CO]  =  nlctrb ( fx, Gx, P) 


O, 

o 

g. 

o 

o, 

o 

o, 

o 

g. 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g. 

o 

o, 

o 

g. 

o 

Q. 

o 

o, 

o 

g, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

g. 

o 

g„ 

o 

g. 

o 

g. 

o 

g„ 

o 

g. 

o 

o, 

o 

g. 

o 

g. 

o 

g, 

o 

g. 

o 


NLCTRB  Compute  the  controllability  matrix  and  associated 

Lie  algebra  for  a  nonlinear  control-affine  system 
(with  or  without  drift)  of  the  following  form: 

xdot  =  f (x)  +  G(x)u,  f(P)  =  0,  x->RNx,  u->RNu 

The  following  three  cases  pertain  for  the  condition  of  the 
dimensions  of  the  Lie  Algebra  Matrix  L  where  Nu  is 
the  dimension  of  the  control  space  and  Nx  is  the 
dimension  of  the  state  space: 

1.  dim (L (Delta) )  =  Nu  The  system  is  completely  integrable 

2.  Nu  <  dim (L (Delta) )  <  Nx  The  system  is  nonholonomic,  but  not 

Small-Time  Locally  Controllable  (STLC) 

3.  dim (L (Delta) )  =  Nx  The  system  is  nonholonomic  and  STLC 

if  the  system  is  driftless  (i.e.  f  0 (x)  =  0  for  all  x  in  X)  enter 
the  zero  vector  for  fO (x) . 

Evaluate  this  script  to  establish  symbolic  objects  for  state 

variables  prior  to  calling  function  where  state  =  number  of 
state  variables 

syms  X 

for  a  =  1: state 

x  =  strcat ( ' x ' , num2str (a) ) ; 
syms (x) ; 

X  ( a )  =  x  ; 

end 


%  This  algorithm  uses  the  P.Hall  Basis  algorithm  as  explained  in 

%  LaValle,  S.M.  Planning  Algorithms,  2006. 

O, 

o 

%Developed  by  Jason  S.  Hall 
%Naval  Postgraduate  School 
%Feb  2007 

%Modified  July  2009  to  work  with  new  Symbolic  Matlab  Toolbox 

S-S-5'9'2'9'9'9'9'9'9'9'5'9'9'9'9'9'5'9'9'9'2'9'2'9'9'9'5'9'9'9'9'9'2'9'5'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'5'9'9'9'5'9'9'9'2'9'9'9'2'9'9'9'5'9' 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%%  State  Variable  Assignment 
[Nx,Nu]  =  size(Gx); 

X  =  sym ( zeros (Nx, 1 )) ; 

for  a  =  1 : Nx 

x  =  strcat (’ x 1 , num2str (a) ) ; 
syms (x); 

X  (a)  =  x ; 

end 


%%  Vector  Variable  Assignment 
syms  h  fO 
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1 


if  isequal ( fx, zeros (Nx, 1 ) )  == 
dimh  =  Nu; 


q  { l } 

=  'hi'; 

V 

=  genvarname ( ' h ' , who) ; 

eval (  [v 

'  =  Gx ( : , 1 )  '  ]  ) 

L  { 1 } 

=  q  { 1 }  ; 

for  a  = 

2 : dimh 

q{a} 

=  strcat (' h ', num2str (a 

V 

=  genvarname ( ' h ' , who) ; 

eval ( [v  '  =  Gx ( : , a) ; ' ] ) ; 

L(a) 

=  q  (a)  ; 

end 

st 

=  1; 

dimh 

=  Nu  +  1; 

q  { l } 

=  '  f  0  '  ; 

f0 

=  fx; 

st 

=  0; 

for  a  = 

2 : dimh 

q{a}  =  strcat ( ' h ' , num2str (a-1 
v  =  genvarname ( ' h 1 , who) ; 

eval ( [v  '  =  Gx ( : , a-1 ) ;  '  ] )  ; 

L (a-1)  =  q (a)  ; 

end 

end 

dim  =  ones ( 1 , Nu+1 ) ; 

Xn  =  zeros (Nx, 1) ; 

for  a  =  1 : Nx 

if  P(a)  ~=  0; 

Xn (a)  =  randn; 

end 

end 

CO  =  Gx; 

Ldel_x  =  double ( subs (Gx, X, Xn) ) ; 
rhmat  =  rank(Ldel  x) ; 

%%  Required  Index  Assignment 
[kk,ll]  =  size(Gx); 

k  =  1; 

j  j  =  i  ; 

%%  Algorithm 
for  a  =  2 : Nu+5 

[o,  p]  =  size (q) ; 
if  rhmat  ==  Nx 
break 

end 

for  b  =  st:Nu 

if  rhmat  ==  Nx 
break 

end 

for  c  =  1 : p 

if  rhmat  ==  Nx 
break 

end 

if  b  ==  0 
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rl  =  ' f 0 ' ; 

else 

rl  =  strcat ( ' h ' , num2str (b) ) ; 

end 

if  iscellstr (q (a-1, c) )  &&  isequal (char (q (a-1, c) ) ,  'n 

r2  =  strcat (q (a-1, c) ) ; 
r  =  strcat ( rl r2 ) ; 
s  =  sym (char (r) ) ; 

else 

r  =  {  '  n '  } ; 
s  =  sym (char (r) ) ; 

end 

if  strcmp ( rl , r2 ) 
r  =  {  '  n  '  }  ; 
s  =  sym (char (r) ) ; 

el seif  (strcmp ( rl ,  ' hi ' )  &&  strcmp ( r2 ,  ' f 0 , hi ' ) )  I  I  . . 

(strcmp (rl, ' h2 ' )  &&  strcmp (r2, ' fO, h2 ') ) 

r  =  {  '  n  '  }  ; 
s  =  sym (char (r) ) ; 

else 

if  size(rl,2)  <=  3 
r3  =  char (r2 ) ; 

if  str2double (rl (2 : end) )  >=  ... 

str2double (r3 (2 : end) )  &&  ... 

isequal (r3 (end) ,' n ' )  ==  0 
r  =  {  '  n  '  }  ; 
s  =  sym (char (r) ) ; 

end 

end 

if  isequal (r, { ' n ' } ) ~=1 
q  (a,  k)  =  r; 
k  =  k  +  1; 

end 

if  length(rl)  >=  2  &&  isequal (char (s) ,' n ' )  ~=  1 

s  2  =  s  ; 

end 

for  d  =  1 : dimh 

if  size(rl)  ==  size (char (q ( 1 , d) ) ) 
if  rl  ==  char(q(l,d)) 
for  e  =  1 : Nu+1 

if  length (s)  ==  e+1 

si  { e }  (dim  (e)  ,  : )  =  s; 
kk  =  dim (e) ; 
dim(e)  =  dim(e)  +  1; 

end 

end 

end 

end 

end 

end 

if  exist  (' s2 ' var ' ) 
fng  =  eval ( s2 ) ; 
zz  =  size(fng,2); 
while  zz  >  1 

zz  =  zz  -1; 


137 


fv  =  fng  (  : , zz) ; 
g  =  fng ( : , zz+1) ; 
dfdx  =  j acobian ( fv, X) ; 
dgdx  =  j acobian (g, X) ; 
fng(:,zz)  =  dgdx*fv-dfdx*g; 

end 

brack  =  fng(:,zz); 

FNG  =  vpa ( subs (brack, X, Xn) , 2 ) ; 

FNG1  =  subs (brack, ' X ' , '  P  '  )  ; 

hmat  =  [Ldel  x,FNG]; 

pp  =  size (hmat, 2 ) ; 

rhmat  =  double (rank (hmat) ) ; 

if  rhmat  ==  pp 

Ldel  x  =  hmat; 


[CO, FNG1 ] ; 

char (si {end} (kk, :)); 
cellstr(qn(9:size(qn,2)-2) ) 
11  +  1; 


CO 

qn 


L (11+1) 
11 


end 

if  rhmat  ==  Nx 
break 

end 

clear  s2 


end 

end 

end 

jj  =  jj+i; 

k  =  1; 

end 

PHall  Basis  =  o; 

if  Nu  <  Nx  &&  isequal ( fx, zeros (Nx, 1 ) )  ==  1 
if  rhmat  ==  Nx 

fprintf('The  system  is  STLC  \n  \n') 
elseif  rhmat  ==  Nu 

fprintf('The  system  is  completely  integrable  \n  \n') 
elseif  rhmat  <  Nx  &&  rhmat  >  Nu 

fprintf('The  system  is  not  STLC  \n  \n') 

else 

fprintf('The  system  doesn''t  fit  \n  \n') 

end 

elseif  Nu  <  Nx  &&  isequal ( fx, zeros (Nx, 1 ) )  ==  0 
if  rhmat  ==  Nx 

fprintf('The  system  is  at  least  accessible,  ') 
fprintf (' verify  no  bad  brackets  to  determine  STLC\n  \n') 
elseif  rhmat  ==  Nu 

fprintf ('The  system  is  completely  integrable  \n  \n') 
elseif  rhmat  <  Nx  &&  rhmat  >  Nu 

fprintf ('The  system  is  not  accessible  \n  \n') 

else 

fprintf ('The  system  doesn''t  fit  \n  \n') 

end 

end 

fprintf (' L (Delta)  =  span  {%s ' , char (L (1) ) ) 
for  a  =  2: length (L) 
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fprintf ( ' , %s ' , char (L (a) ) ) 

end 

fprintf ('}  \n ' ) 

fprintf (' \nThe  LRAC  is  %i  \n',  rhmat) ; 

fprintf (' P . Hall  Basis  to  Depth  %i  \n',PHall  Basis) 
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APPENDIX  B.  PROOFS  OF  EMPLOYED  LEMMAS 


Lemma  1. 

Proof. 


to 


(?)  is  uniformly  bounded,  i.e.  to a(t)  <c\/t 


1  2 

By  rearranging  V  =  —  k  ||coa  ||  +  2  (l  -  q4 )  to  yield 


^|to„||2=2L-4(l-^4) 


(255) 


and  then  substituting  this  into  V  =  qu(j){xu )  - aq]  - kd  ||co0 1|2  produces 

v  =  -2  dV  +  4d  (1  -  q.)  +  )  -  aq:  (256) 

The  desired  stable  linear  dynamics  ^(xi()  +  2y0(xH)  +  y2^(xi()  =  0  and  the 
condition  det[yo(x)]*0  for  the  infinite  set  of  control  laws  +Pd  (x)y  with 

=  a^  (x)h(x)  ensure  that 

lim^(x  )  =  0  (257) 

f— >oo  x  7 


for  an  arbitrarily  small  selection  of  that  damping  coefficient  /)  on  the  controls  coefficient 
generalized  inverse  a)(x).  Therefore,  given  c/4|  <  1  Vt  from  the  definition  of  the  error 
quaternion  (Cristi  and  Burl  1993),  F(t)  is  bounded  and  ||<oa  (t)||  is  bounded. 

Lemma  2.  If  0(xM)(t)  — >0  exponentially,  i.e.  )(t)|  <  Ae~rt ,  t  >  0  for  some  A,  y  >  0 , 
then 

lim/(t)  =  /(oo)<+oo  (258) 

f-H-  oo  V  7  V  7 

where  f  (  t)  is  defined  as 
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(259) 


/W=jK(T)|2rfr 

0 


Proof.  For  all  fel  (Cristi  and  Burl  1993) 


V(t)  =  V(ti)  +  \v{r)dT  (260) 

0 


1  ||2 

Substitution  of  the  Lyapunov  function  V  given  by  V  =  —k  |(Oa  |  +  2  (l  -  q4  )  and  its 
derivative  V  =  qu</>{*u)- ccq 2  - kd ||<oa ||2  yields 


v(t)  =  v(0i)-kd\\wa(T)\  dT  +  \qXe(T)(j)(xu)(T)dT-a\ql(T)dT  (261) 


Since  F(7)  >  0  Vt 


t  t 

kd\\ (Oa  (r)||“  dr  <  V (0)  +  J|^(xH  (262) 

0  0 


Letting  A  =  Xjkd  yields 


JIK  M||2  dr<AV(ti)  +  a||^(xb)(t)|  dr 
0  0 


(263) 


and  since  ^(xu)  exponentially  decays  through  the  selection  of  y  in 
(/)(\u^  +  2y(/>  ( \u)  +  y2(/>  ( xH )  =  0  to  satisfy  the  characteristic  equation  given 
by  s2  +  2 ys  +  y2  =  0  ,  the  integral  is  finite.  Then  for  some  constant  c 

t 

f(t)  =  J|(Oa  (r)||~  dr  <  c  Vt  (264) 

0 

Given  (264)  is  monotonically  non-decreasing  and  upper  bounded,  it  has  a  finite 
limit  as  t  +oo  and  thus  for  some  constant  c 
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-I-UU 

fH=  JKWlf dT- 


(265) 


Lemma  3.  Given  (257), 


limU{xu)\\  =  a>a+aqu=0 

t  — ^co  v  7 


(266) 


which  implies  that 


\vm\(Ou{t)\  =  -aqu{t) 


(267) 


Furthermore,  with  the  lim  o>  (?)  =0  from  (250)  and  the  kinematic  equation  for 

t—>oo  \ 


the  scalar  component  of  the  error  quaternion  given  by  qA-  0  xu  +  Hp^a 


*-*  ®  2 


(268) 


/'.y  implied.  Therefore,  substitution  of  (267)  into  (268)  yields 


lim  1^4 1  =  -“C 

f->co  2 


(269) 


so  that 


lim  q4  (?)  =  0 

oo  V  7 


(270) 


Additionally,  given  the  kinematic  equation  for  the  unactuated  component  of  the 
error  quaternion  expressed  by  qu  =  ^{q4mu  -  q^G>fl  )  ,the  lim||(Da  (?)||  =  0  from  (250)  and 
the  result  from  (270)  yields 


lim  qu  (?)  =  0 


(271) 


And,  therefore,  from  (267), 


143 


(272) 


lim \co  (?)  =0 

and  furthermore 

limllx  (?)||  =  0  (273) 

?— >oo  >1  V  7  II 

Lemma  4.  Given  lim  1 1©  (?)||  =  0  from  (250)  and  limllx  (?)||  =  0  from  (273),  and  given 
the  kinematic  equations  for  the  error  quaternion  as  expressed  by 

*1  p~  ®2xl  ^21  )  X«  ^22 

lim|k(?)|  =  0  (274) 

Proof.  From  qp=[o2ll  F21  (q^)]x„  +F22  (q*^)©,,  given  lim||©a  (?)||  =  0 

from  (250)  and  limllx  (?)||  =  0  from  (273) 

t — >oo  V  7  II 

HM?)l=o  ^275) 

t— >00  II  y  II 

(oa  =  -<7(oa  -  A'q/;  can  be  rewritten  as 

(®« +  °<\P ) ++„  +  P )  =  °<\P  (276) 

where  k  =  d  >  0,  cr  =  kd  l  >  0 .  By  defining 

e  (?)  =  ©a  (?)  +  crq^  (? )  (277) 

and  substituting  this  into  (276)  yields 

e(?)  =  -A-e(?)  +  crqp(?)  (278) 

From  (275)  and  given  k  >  0 

lim|e(?)||  =  0  (279) 

t — >oo  H  V  7H 
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and,  therefore,  substitution  of  (279)  into  (277)  yields 


lim  q  (f)  =0 


(280) 
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APPENDIX  C.  MATLAB  CODE  FOR  UNDERACTUATED 
QUATERNION  FEEDBACK  CONTROLLER 


function  [T]  =  RateStable  (w,  q,  qd,  wd,  J,  d,  k,  gamma,  ... 

Betal,  Beta2,  al,  select) 

%  Quaternion  Feedback  Attitude  Regulator  for  Underactuated  Spacecraft 


for 


end 

for 


i  =  2:4 
if  q ( i )  == 
q  ( i )  = 

end 


0 

Betal ; 


i  =  1:3 
if  w(i)  == 


==  0; 


w(i)  =  Betal; 

end 

end 

q0 

= 

q  ( 1 )  ; 

qi 

= 

q  ( 2 )  ; 

q2 

= 

q  ( 3 )  ; 

q3 

= 

q  (4)  ; 

qdO 

= 

qd ( 1 ) ; 

qdl 

= 

qd ( 2 ) ; 

qd2 

- 

qd ( 3 ) ; 

qd3 

= 

qd ( 4 ) ; 

qdv 

= 

[qdl ; qd2 ; qd3 ] ; 

qv 

= 

[ql;q2;q3] ; 

qx 

= 

[0  -q3  q2;q3  0  -ql;-q2  ql 

0]  ; 

qOe 

= 

q0*qd0  +  qv'*qdv; 

qve 

= 

qd0*qv  -  q0*qdv  +  qx*qdv; 

qle 

= 

qve ( 1 ) ; 

q2e 

= 

qve ( 2 ) ; 

q3e 

= 

qve ( 3 ) ; 

we 

= 

w  -  wd ; 

wle 

= 

we  ( 1 )  ; 

w2e 

= 

vie  (2)  ; 

w3e 

= 

vie  (3)  ; 

w  X 

= 

[0  -w3e  w2e;w3e  0  -wle;-w2e  wle 

0 

Jll 

= 

J(l, 1) ; J12  =  J ( 1 , 2 ) ; J13  = 

J  ( 1 , 3 

)  ; 

J21 

= 

J (2 , 1) ; J22  =  J (2 , 2 ) ; J23  = 

J  (2 , 3 

)  ; 

J31 

= 

J  (3, 1)  ;  J32  =  J  ( 3 , 2 )  ;  J33  = 

J  (3, 3 

)  ; 

a 

_ 

[1/2* (-2*wle*Jll*J13*J22+2 

:*wle*. 

Ji: 

2*wle*J31*J22*J33-2*J12*w3e*J13*J32-2*J22*w3e*J23*J32- . 
2*J22*wle*J12*J23-2*J32*wle*J12*J33+2*w3e*J13*J12*J23+. 
2*w3e*J33*J23*J32+2*w3e*J12A2*J33+2*w3e*J22A2*J33+. . . 
2*wle*J13*J22A2+2*wle*J13*J32A2-2*w3e*J13A2*J22-. . . 
2*w3e*J22*J33A2+4*w2e*J12A2*J23+4*w2e*J23*J32A2- . . . 
4*w2e*J12*J13*J22-4*w2e*J32*J22*J33+. . . 
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Lf 2phi 


2* 

wle* 

J31* 

J23* 

J32- 

al*q3e* J1 1 

al 

*q3e 

*  J2 1 

*  J12 

*  J33 

-al*q3e* J2 

al 

*q3e 

*  J31 

*  J13 

*  J22 

) / ( Jll* J22 

J2 

1* J13* J32+J3 

1* J12* J23- J31* 

1/ 

2*  (2 

*wle 

*J11 

*  J12 

* J33-2  *wle 

2* 

wle* 

J2 1  * 

J22* 

J33- 

2*wle*J21* 

2* 

w2e* 

J22  A 

2  * J33-2  * 

w2e* J13 A2* 

2* 

w2e* 

J22  * 

J33  A 

2-4* 

w3e* J13 A2* 

2* 

w2e* 

J12  * 

J13* 

J32- 

2*w2e* J22* 

2* 

J23* 

wle* 

J13* 

J22+2* J33*w2e* 

2* 

J33* 

wle* 

J13* 

J32+4*w3e*J13* 

al 

*q2e 

*J11 

*  J22 

*  J33 

-al*q2e* J1 

al 

*q2e 

*  J2 1 

*  J13 

* J32+al*q2e* J3 

al 

*q2e 

*  J31 

*  J13 

*  J22 

) / ( J11*J22 

J2 

1* J12* J33+J2 

1* J13* J32+J31* 

1/2 

*al* 

wle* 

(-1/ 

2  *ql 

e*wle-l/2* 

l/2*al*w3e* (1 
l/2*al*w2e*  (- 
(l/2*al 
J11*J23 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
J1 
J3 
(J1 
J31 
(J1 
J31 
(J1 
J31 
(J1 
J31 
(J1 
J31 
(J1 
J31 
(J1 
J31 


*J22*J33+al*q3e*Jll*J23*J32+. 
l*J13*J32-al*q3e*J31*J12*J23+ 
*J33-J11*J23*J32-J21*J12*J33+ 
J13* J22 ) ; 

* J1 1  * J13* J32+ . . . 
J23*J32+2*w2e*J12A2*J33+. . . 
J22-2*wle*J12*J23A2- . . . 
J32-4*w3e*J23A2*J32- . . . 
J23*J32+2*J13*w2e*J12*J23+. . . 
J23*J32-2*wle*J12*J33A2+. . . 
J12*J33  +  4*w3e* J23*J22*J33+ . . . 
l*J23*J32-al*q2e*J21*J12*J33+ 
1* J12* J23- . . . 
*J33-J11*J23*J32- . . . 
J12*J23-J31*J13*J22) ] ; 
q2e*w2e-l/2*q3e*w3e)+. . . 
*q0e*w2e-l/2*qle*w3e) - . . . 
2*qle*w2e+l/2*q0e*w3e)+. . . 


/2*q3e*wle+l/2 
l/2*q2e*wle+l/ 

*qOe+ (- (J12*J23-J13*J22) / ( J11*J22* J33- . . . 
*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- .  .  . 
1*J13*J22) * J21- (J12*J33-J13*J32) / ( J11*J22*J33- . . . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- .  .  . 

1* J13* J22 ) * J31 ) *wle+ ( ( J12 * J33- J13* J32 ) / ( Jll* J22* J33- . . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 

1* J13* J22 ) *w3e+ ( J12* J23- J13* J22 ) / ( Jll* J22* J33- . . . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 

1* J13* J22 ) *w2e) *J11+ ( ( J22*J33-J23*J32) / ( Jll* J22* J33- . . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 
1*J13*J22) *w3e- ( J12* J2  3- J13* J22 ) / ( J11*J22*J33-  .  . . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 

1* J13* J22 ) *wle) * J21+ (- ( J22* J33- J23* J32 ) / ( Jll* J22*J33- . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- .  .  . 

1* J13* J22 ) *w2e- ( J12 * J33- J13* J32 ) / ( Jll* J22* J33-  .  .  . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 

1* J13* J22 ) *wle) * J31+ (- ( J12*J23-J13*J22) / ( Jll* J22* J33- . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23-. . . 
1*J13*J22) *J22- ( J12  * J33- J13* J32 ) / ( Jll* J22* J33- .  .  . 
1*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 
1*J13*J22) * J32 ) *w2e+ ( - ( J12* J23- J13* J22 ) / ( Jll* J22*J33- . 
1*J23*J32-J21*J12*J33+J21*J13*J32+.  .  . 
1*J12*J23-J31*J13*J22) *J23- ( J12* J33- J13* J32 ) /. . . 
1*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 

* J12* J23- J31* J13* J22 ) *J33) *w3e) * ( ( ( ( J12* J33- J13* J32 ) / 
1*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
*J12*J23-J31*J13*J22) *w3e+ ( J12* J23- J13* J22 )  /.  .  . 
1*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
*J12*J23-J31*J13*J22) *w2e) *J11+ ( ( J22*J33-J23*J32)  /.  .  . 
1*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+.  .  . 
*J12*J2  3-J31*J13*J22) *w3e- ( J12* J23- J13* J22 )  /.  .  . 
1*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+. .  . 
*J12*J23-J31*J13*J22) *wle) *J21+ (- ( J22* J33- J23* J32 ) / . . 
1*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
*J12*J23-J31*J13*J22) *w2e- ( J12* J33- J13* J32 ) /. . . 
1*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
*J12*J23-J31*J13*J22) *wle) *J31) *wle+ . . . 
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( ( (J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32- .  .  . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. . . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J12+. . . 

( (J22*J33-J23*J32)  / ( Jll* J22* J33-J11* J23* J32-J21* J12* J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- . . . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J22+. . . 

(- ( J22* J33- J2  3* J32 ) / ( J11*J22*J33-J11* J23* J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e- . . . 
(J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J32) *w2e+ . . . 

( ( (J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+ . . . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J13+. . . 

( (J22*J33-J23*J32) / ( Jll* J22* J33-J11* J23* J32-J21* J12* J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- .  .  . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J23+. . . 

(- ( J22* J33- J2  3* J32 ) / ( J11*J22*J33-J11* J23*J32-  .  .  . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e- . . . 
(J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J33) *w3e) +. . . 
(-l/2*al*q3e+ ( ( J12*J23-J13*J22) / ( Jll* J22* J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 
J31*J13*J22) *J11- (J22*J33-J23*J32) / ( Jll* J22* J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- .  .  . 
J31*J13*J22) * J31 ) *wle+J23-J13*J22) / ( J11*J22*J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 
J31*J13*J22) * J12- ( J22* J33- J23* J32 ) /  ... 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *J32) *w2e+ ( ( J12 * J33- J13* J32 )  /  .  .  . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *w3e+ ( J12* J23- J13* J22 )  /.  .  . 
(J11*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+.  .  . 
J31*J12*J23-J31*J13*J22) *w2e) *J12+ ( ( J22*J33-J23*J32) /. . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *w3e- ( J12* J23- J13* J22 ) /. . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *wle) *J22+ (- ( J22* J33- J23* J32 ) / . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J2  3-J31*J13*J22) *w2e- ( J12* J33- J13* J32 )  /.  .  . 
(J11*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+. .  . 
J31*J12*J23-J31*J13*J22) *wle) *J32+ ( ( J12*J23-J13*J22) /. . . 
(J11*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+. .  . 
J31*J12*J23-J31*J13*J22) *J13- ( J22* J33- J23* J32 )  /.  .  . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *J33) *w3e) * ( ( (- (J11*J33- . . . 

J13* J31 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- (J11*J23- . . . 
J13*J21) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J11+. . . 

(- ( J21*J33-J23*J31) / (J11*J22*J33-J11*J23*J32- . . . 
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J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. . . 
(J11*J23-J13*J21) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J21+. . . 

( (J21*J33-J23*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e+ (J11*J33- . . . 

J13* J31 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+.  .  . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J31) *wle+ . . . 

( (- ( Jll* J33- J13* J31 ) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- . . . 
(J11*J23-J13*J21) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J12+. . . 

(- (J21*J33-J23*J31) / (J11*J22*J33-J11*J23*J32-  .  .  . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. . . 
(J11*J23-J13*J21) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+.  . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J22+. . . 

( ( J21* J33- J23* J31 ) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e+. . . 
(J11*J33-J13*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J32) *w2e+. . . 

( (- (J11*J33-J13*J31) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- . . . 
(J11*J23-J13*J21) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J13+. . . 

(- ( J21* J33- J23* J31 ) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. . . 
(J11*J23-J13*J21) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J23+. . . 

( (J21*J33-J23*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e+. . . 
(J11*J33-J13*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J33) *w3e) +. . . 
(l/2*al*q2e+ ( ( J12 * J33- J13* J32 ) / (J11*J22*J33-J11*J23*J32- 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *J11+. . . 

( J22* J33- J23* J32 ) / ( J11*J22*J33-J11*J23*J32-J21* J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *J21) *wle+ . . . 

( (J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *J12+. . . 

( J22* J33- J23* J32 ) / ( J11*J22*J33-J11* J23*J32-J21* J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *J22) *w2e+. . . 

( (J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J2  3-J31*J13*J22)  *J13+.  .  . 

( J22* J33- J23* J32 ) / ( Jll* J22* J33-J11* J23* J32-J21* J12* J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *J23) *w3e+ . . . 

( (J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. . . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J13+. . . 

( (J22*J33-J23*J32) / ( Jll* J22*J33-J11*J23*J32-J21*J12* J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- . . . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J23+. . . 

(- (J22*J33-J23*J32) / ( Jll* J22* J33-J11* J23* J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e- . . . 
(J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J33) * . . . 
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( ( ( (J11*J32-J12*J31) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+. . . 
(J11*J22-J12*J21) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J11+ ( (J21*J32- 
J22*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- (J11*J22- . . . 

J12 * J2 1 )  / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J21+. . . 

(- ( J21* J32- J22* J31 ) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e- . . . 
(J11*J32-J12*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13* J32+J31* J12* J23- . . . 

J31* J13* J22 ) *wle) *J31) *wle+ ( ( ( J11*J32-J12*J31) / . . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *w3e+ ( J11*J22-J12*J21)  /  .  .  . 
(J11*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+.  .  . 
J31*J12*J23-J31*J13*J22) *w2e) *J12+ ( ( J21*J32-J22*J31)  /.  .  . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J2  3-J31*J13*J22) *w3e- ( J11*J22-J12*J21)  /.  .  . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *wle) *J22+ (- ( J21*J32-J22*J31) / . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *w2e- ( J11*J32-J12*J31) /. . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *wle) *J32) *w2e+ ( ( (J11*J32- . . . 
J12*J31 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+ (J11*J22- . . . 

J12  * J2 1 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J13+. . . 

( (J21*J32-J22*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+ 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- (J11*J22- . . . 

J12  * J2 1 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J23+. . . 

(- ( J21* J32- J22* J31 ) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e- . . . 
(J11*J32-J12*J31) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J33) *w3e) ; 
Lfphi  =  al* (l/2*q0e*wle-l/2*q3e*w2e+l/2*q2e*w3e) + ( ( (J12*J33- . . . 

J13* J32 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+ . . . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23- .  .  . 

J31* J13* J22 ) *w2e) *J11+ ( ( J22*J33-J23*J32) / . . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *w3e- ( J12*J23-J13*J22) /. . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *wle) *J21+ (- ( J22*J33-J23*J32) / . . 
(J11*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+.  .  . 
J31*J12*J2  3-J31*J13*J22) *w2e- ( J12* J33- J13* J32 )  /.  .  . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J2  3-J31*J13*J22) *wle) *J31) *wle+ ( ( (J12*J33- .  .  . 
J13* J32 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e+ . . . 
(J12*J23-J13*J22) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. 
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J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e) *J12+. .  . 

( (J22*J33-J23*J32) / ( Jll* J22* J33-J11* J23* J32-J21* J12* J33+ . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- (J12*J23-  .  .  . 

J13* J22 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J22+.  .  . 

(- ( J22* J33- J23* J32 ) / ( Jll* J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23-J31*J13*J22) *w2e-. . . 
(J12*J33-J13*J32) / (J11*J22*J33-J11*J23*J32-  .  .  . 
J21*J12*J33+J21*J13*J32+J31*J12*J2  3- . .  . 

J31*J13*J22) *wle) *J32) *w2e+ ( ( ( J12 * J33- J13* J32 ) / . . . 
(J11*J22*J33-J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *w3e+ ( J12*J23-J13*J22)  /  .  .  . 
(J11*J22*J33-J11*J2  3*J32-J21*J12*J33+J21*J13*J32+.  .  . 
J31*J12*J2  3-J31*J13*J22) *w2e) *J13+ ( (J22*J33- .  .  . 

J23* J32 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. .  . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *w3e- (J12*J23- .  .  . 

J13* J22 ) / (J11*J22*J33-J11*J23*J32-J21*J12*J33+. . . 
J21*J13*J32+J31*J12*J23-J31*J13*J22) *wle) *J23+. . . 

(- ( J22* J33- J23* J32 ) / ( Jll* J22*J33-J11*J23*J32- . . . 
J21*J12*J33+J21*J13*J32+J31*J12*J23- . . . 

J31* J13* J22 ) *w2e- ( J12* J33- J13* J32 ) / ( J11*J22*J33- . . . 
J11*J23*J32-J21*J12*J33+J21*J13*J32+. . . 
J31*J12*J23-J31*J13*J22) *wle) *J33) *w3e; 
phi  =  wle  +  al*qle; 

phidot  =  Lfphi; 

b  =  -Lf2phi  -  2 *gamma*Lfphi  -  gammaA2*phi; 

wu  =  wle; 

wa  =  [w2e;w3e] ; 

qa  =  qve (2:3) ; 

%qu  =  qve  ( 1 )  ; 

Fw  =  -inv ( J) *wx* J; 

521  =  Fw (2:3,1) ; 

522  =  Fw  (2  :  3, 2  :  3)  ; 

if  sqrt(a'*a)  >=  Betal 
adinv  =  a/(a'*a); 
else  adinv  =  a/BetalA2; 
end 

Pd  =  eye (2)  -  adinv*a'; 

udbar  =  adinv*b; 
if  select  ==  1 

%%  Lyapunov  Function-Based  Control  Law 
c  =  wa'*Pd*wa; 
if  sqrt(c)  >=  Beta2 

etad  =  wa/ (wa ' *Pd*wa) ; 
else  etad  =  wa/Beta2; 
end 

y  =  etad* (-wa ' *S21*wu  -  wa'*S22*wa  -  wa'*udbar  -  k*qve'*we)  -  d*wa 

else 

%%  Feedback  Linearizing  Control  Law 
y  =  (-S21*wu  -  S22*wa  -  k*qa  -  d*wa) ; 

end 
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tau_a 

tau 

T 


=  udbar  +  Pd*y; 
=  [0;tau_a]; 

=  J*tau; 
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